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Demand  for  advances  in  the  theoretical  and  experimental  understanding  of 
high  performance  materials  has  provided  the  impetus  to  organize  this 
Workshop  on  ‘Composite  Material  Response:  Constitutive  Relations  and 
Damage  Mechanisms’.Nlt  was  held  at  the  Stakis  Grosvenor  Hotel  in 
Glasgow  from  July  30  w  31,  1987.  Experts  from  six  different  countries 
including  Australia,  Greece,  Italy,  Poland,  United  Kingdom  and  United 
States  were  invited  to  present  their  most  recent  research  findings  and  ideas 
on  composite  material  bdiiaviour.  Particular  emphases  were  placed  on  the 
damage  mechanisms  associated  with  mechanical,  thermal  and/or  chemical 
effects.  Their  influence  on  the  ways  with  which  constitutive  relations 
are  formulated  is  relevant  in  quantifying  the  behaviour  of  application- 
specific  materials  such  as  composites.  Such  a  knowledge  is  particularly 
lacking. 

The  need  for  more  effective  use  of  advanced  materials  was  emphasized  in 
the  Opening  Address  by  Dr  Fritz  H.  Oertel,  Jr,  who  is  Chief  of  Research 
in  the  US  Army  Research  Development  and  Standardization  Group — 
United  Kingdom  (USARDSG-UK).  This  Group  grew  out  of  the  American, 
British,  Canadian  and  Australian  (ABCA)  war-time  alliance  and  sponsors 
research  activities  to  the  interests  of  the  US  Army  and  the  appropriate 
European/Middle  Eastern/ African  Army  or  Institution.  To  be  acknow- 
I  ledged,  in  particular,  is  the  support  of  this  Workshop  provided  by  the 
^USARDSG-UK  Physics/Mathematics  Branch  headed  by  Dr  Julian  J.  Wu. 
It  provided  an  opportunity  for  evaluating  the  current  as  well  as  future  needs 


in  advanced  materials  science  and  technology. 
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Microstructure  and  Damage  Dependence  of 
Advanced  Composite  Material  Behavior 

G.  C.  SiH 

Institute  of  Fracture  and  Solid  Machanics,  Lehigh  University, 
Bethlehem,  Pennsylvania,  USA 


ABSTRACT 

The  macromechanical  behavior  of  composite  materials  exhibits  complex 
dependence  not  only  on  the  microstructure  but  also  in  the  way  with  which 
damage  occurs.  This  is  a  load-time  history  dependent  process,  regardless  of 
whether  composites  are  being  used  for  their  many  unique  behaviors  in  terms  of 
strength,  fracture  resistance,  thermal  properties,  etc.  Because  the  prevention 
of  mechanical  failure  is  almost  always,  an  important,  if  not  critical, 
requirement,  damage  accumulated from  delamination,  initiation  and  growth  of 
cracks  in  the  fiber  and  matrix  and  at  the  interface  becomes  an  important 
consideration  as  the  sequence  of  failure  modes  can  affect  the  outcome.  Many  of 
the  material  testing  procedures  developed  for  metals  may  not  be  valid  because 
each  specimen,  even  though  loaded  uniaxially,  possesses  a  unique  behavior  on 
account  of  its  highly  nonhomogeneous  microstructure  and  damage  pattern. 
Usage-specific  is,  in  fact,  a  salient  feature  of  composites  that  can  be  tailor- 
made  to  satisfy  certain  performance  requirements. 

The  prediction  andjor  quantitative  assessment  of  composite  material 
behavior  remain  undeveloped.  This  is  not  surprising  because  the  classical 
theories  of  mechanics  and  physics  are  not  conducive  for  analysing  inhomo¬ 
geneity  arising  from  the  material  microstructure  and  nonuniform  damage. 
Notwithstanding  such  complexities,  this  communication  will  attempt  to 
address  a  basic  feature  of  composite  behavior,  i.e.  nonuniform  energy  dis¬ 
sipation  as  a  result  of  local  temperature  fluctuation.  This  is  also  related  to 
change  in  local  strain  rates  and  strain  rate  history  being  the  direct  cause  of 
microstructure  inhomogeneity  and  damage.  Load  transfer  characteristics  as 
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altered  by  processing  across  fiber/matrix  interfaces  or  bonding  surfaces  in 
laminates  can  alter  the  composite  properties  even  when  the  same  fibers  and 
matrix  materials  are  used.  There  is  the  exchange  of  surface  and  volume  energy 
across  an  interface  that  involves  the  interaction  of  mechanical  and  physico¬ 
chemical  effects.  The  fluctuation  of  the  local  volume  energy  density  dW/dV 
and  surface  energy  density  fTNjdK  in  composites  would,  therefore,  be 
regarded  as  the  quantities  of  interest.  They  are  related  by  the  rate  change  of 
volume  and  surface  dW/d\  that  serves  as  a  scaling  length  parameter  and  an 
index  related  to  the  degree  of  homogeneity  or  inhomogeneity  of  the  system. 
Here,  homogeneity  includes  the  combined  effect  of  material  microstructure, 
geometry  and  load  type.  Illustrative  examples  will  be  provided  and  discussed 
briefly  in  connection  with  the  technique  of  Electromagnetic  Discharge 
Imaging  (EDI ),  a  possible  means  for  evaluating  the  mechanical,  thermal  and 
chemical  effects  associated  with  composite  material  behavior. 


1.  INTRODUCTION 

As  modern  structures  and  vehicles  are  required  to  perform  under  con¬ 
ditions  of  high  strength  and  toughness,  combined  with  low  weight, 
composites  are  becoming  the  material  of  choice.  By  dispersing  particles  or 
fibers  of  one  substance  in  a  matrix,  or  binder,  of  another,  the  constituent 
elements  or  microstructures  can  be  tailor-made  to  meet  the  performance 
requirements.  More  and  more  of  the  components  in  military  and  commer¬ 
cial  aircraft,  automobiles  and  sports  equipment  are  made  of  composites. 
Resistance  to  high  temperature  is  another  added  feature  of  composites  that 
is  so  important  to  the  design  of  rocket-motor  components  and  missile  nose 
cones. 

Composite  behavior  is  complex  because  it  depends  on  how  the  individual 
constituents  are  combined.  The  overall  mechanical/thermal/electrical 
properties  can  differ  widely  even  when  composites  are  made  from  identical 
raw  materials  but  processed  differently.  This  is  because  processing  directly 
affects  the  structure  of  the  final  product.  In  this  context,  structure  refers  to 
geometric  discontinuities  such  as  defects  and  voids  and  molecular  features. 
Polymer  matrix,  for  example,  consists  of  small  molecules  joined  chemically 
in  chains  or  networks  of  repeating  units  to  form  huge  molecules  with  an 
infinite  variety  of  possible  three-dimensional  structures.  The  size,  internal 
arrangement,  chemical  connectivity  and  spatial  distribution  of  the  mole¬ 
cules  are  controlled  by  processing.  It  is  vital,  therefore,  to  be  able  to  quantify 
the  sequence  of  events  which  culminate  in  composite  material  damage.  This 
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requires  a  generalized  approach  that  can  uniquely  and  consistently  identify 
results  at  the  atomic,  microscopic  and  macroscopic  scale.  Material  testing 
methods  developed  for  metals  that  invoke  homogeneity  and  isotropy  are 
not  adequate.  Because  a  major  portion  of  the  useful  life  of  composites 
involves  subcritical  damage  in  the  form  of  fiber  breaking,  matrix  cracking, 
interface  debonding,  delamination,  etc.,  these  mechanisms  of  failure  cannot 
be  ignored  in  predicting  performance.  Composite  specimens  will  respond 
differently,  depending  on  the  combined  interactions  of  load,  size  and 
geometry,  in  addition  to  the  inhomogeneous  nature  of  the  microstructure. 

All  substances  are,  in  fact,  inhomogeneous  and  anisotropic  when  they  are 
examined  at  the  microscopic  level.  It  is  the  sensitivity  of  their  microstruc¬ 
ture  in  responding  to  load  that  determines  whether  details  at  the  lower  scale 
level  need  to  be  analysed  or  not.  While  homogeneity  can  be  assumed  for 
most  polycrystals  deformed  under  normal  loading  conditions,  the  details  of 
their  microstructures  become  increasingly  more  important  when  the 
loading  rate  is  slowed  down  to  that  of  creep.  (The  grains  in  a  structural  steel 
are  made  sufficiently  small  in  comparison  with  the  specimen  size  so  that  the 
influence  of  microstructure  entities  can  be  adequately  reflected  via  the 
macroparameters  such  as  yield  strength,  ultimate  strength,  etc.  The  same 
can  be  done  for  composites  if  the  size/time  influence  of  the  constituents  are 
appropriately  suppressed  for  the  given  loading  rates  under  consideration.) 
Polycrystal  metal  structures  are  no  less  complicated  than  those  in 
composites.  Precise  characterization  of  present-day  composites  that  are 
designed  to  aiter  their  global  response  by  microstructurc  effects  necessitates 
a  knowledge  of  the  way  with  which  energy  dissipates  throughout  the 
inhomogeneous  structure  as  it  is  being  physically  damaged.  This  involves 
the  continuous  fluctuation  of  temperature  as  a  result  of  the  constant  change 
of  volume  with  surface  area,  dV/dA,  for  each  composite  element,  an  effect 
that  has  been  neglected  in  the  development  of  continuum  mechanics 
theories  such  as  elasticity,  plasticity,  etc.  (The  assumption  of  letting  dV/dA 
vanish  in  the  limit  has  been  invoked  in  previous  works’  '  ’’  dealing  with  the 
failure  of  composite  systems.  This  has,  in  fact,  excluded  the  thermal/ 
mechanical  interaction  or  the  mechanism  of  energy  dissipation.  Conven¬ 
tional  continuum  mechanics  theories  are  unable  to  predict  the  pheno¬ 
menon  of  cooling/heating* '  “  in  specimens  or  structures  as  the  mechanical 
load  is  increased  monotonically.  Classical  thermodynamics  are  also  not 
valid  because  reversal  of  heat  transfer  leads  to  a  sign  change  in  entropy.) 

In  order  to  distinguish  the  difference  between  the  rates  at  which  energy  is 
dissipated,  say  in  a  fiber  and  matrix,  it  is  necessary  to  consider  the  exchange 
of  surface  and  volume  energy.  Thermal/mechanical  interactions  across  the 
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interface  or  near  a  defect  also  give  rise  to  temperature  change  and  high  rates 
of  energy  dissipation.  To  be  emphasized  in  this  work  is  the  consumption  of 
energy  in  the  composite  damage  process  that  is  not  only  time  dependent  but 
also  varies  from  location  to  location.  The  size/time/temperature  effect  can 
obscure  the  observation  of  composite  material  behavior  when  the  same 
event  is  examined  at  the  atomic,  microscopic  and  macroscopic  scale.  What 
appears  to  be  cooling  and  dilatation  on  one  scale  level  may  correspond  to 
heating  and  distortion  on  another  level.  Unless  the  fundamental  hurdle  of 
material  characterization  is  overcome,  no  confidence  can  be  placed  in 
predicting  the  performance  of  composites.  Composite  damage  analysis  also 
involves  the  nondestructive  evaluation  of  thermal,  mechanical  and 
chemical  effects.  To  this  end,  advancements  have  been  made  on  the 
Electromagnetic  Discharge  Image  (EDI)  technique‘s-  using  high  voltage 
electrical  discharge.  This  method  is  particularly  suited  for  analysing  failure 
in  composites.  Changes  in  the  microstructure,  chemical  composition,  and 
interface  condition  can  be  identified  with  the  mechanical  behavior  of 
composites  so  that  a  more  precise  distinction  can  be  made  between  the 
intrinsic  and  apparent  material  parameters  at  a  given  scale  level.  A  primary 
function  of  the  EDI  method  is  the  evaluation  of  energy  dissipated  in  various 
different  forms. 


2.  THERMAL/MECHANICAL  INTERACTION: 

FIBER,  MATRIX  AND  INTERFACE 

The  thermal/mechanical/electrical  behavior  of  composites  is  very  much 
a  function  of  the  reinforcing  materials  and  a  synergy  between,  say,  the  fibers 
and  matrix.  When  the  fibers  are  stretched  in  the  absence  of  the  matrix,  they 
will  react  individually.  That  is,  the  fibers  once  broken,  are  unable  to  support 
the  load.  If  the  fibers  are  embedded  in  a  matrix  that  is  more  readily 
deformable,  load  transfer  from  the  fibers  to  the  matrix  can  prevail  even 
when  the  fibers  are  broken.  The  choice  of  a  matrix,  therefore,  determines 
not  only  how  a  composite  must  be  fabricated.  An  understanding  of  the  load 
transfer  characteristics  of  the  individual  components  in  a  composite  is 
necessary  before  examining  the  behavior  of  the  complex  system. 

2.1.  Fiber 

There  are  presently  many  materials  that  can  be  made  in  the  form  of  fibers. 
Elements  such  as  carbon,  aluminum,  silicon,  etc.,  can  be  used  to  form 
compounds  in  which  atoms  are  joined  by  strong  and  stable  bonds.  These 
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Fig.  1.  Cooling/heating  of  606I-T6  aluminum  specimen  in  tension  at  a  displacement  rate  of 

8-467  X  10”’  m/s. 


compounds  will  exhibit  different  thermal/mechanical  properties  when  they 
are  loaded  mechanically.  Figure  1  shows  the  temperature  fluctuation,  i.e. 
0  as  a  function  of  time  t  for  a  6061-T6  aluminum  specimen  extended  with 
a  displacement  rate  of  8-467  x  10'*  m/s.  The  specimen  gauge  length  is 
approximately  57  mm  and  the  cross-sectional  area  is  3-81  mm^.  Contrary  to 
the  ordinary  notion  that  the  material  would  heat  up  when  loaded,  it 
actually  cools  before  returning  to  the  ambient  condition.  The  recovery  time 
is  approximately  26  s  and  is  loading  rate  dependent.  This  means  that  energy 
dissipation  is  strictly  rate  dependent.  It  has  been  shown  in  Ref.  1 1  that  the 
onset  of  heating  can  be  delayed  to  200  s  if  the  loading  rate  is  reduced  by  one 
order  of  magnitude.  Dissipation  has  been  known*  “  “  to  be  small  during 
cooling  and  it  increases  at  an  extremely  high  rate  when  heating  starts. 
Therefore,  mismatch  of  energy  dissipation  rates  between  the  fiber  and 
matrix  can  be  undesirable  as  it  tends  to  decrease  the  effectiveness  of  load 
transfer. 

2.2.  Matrix 

If  the  same  aluminum  specimen  discussed  earlier  was  made  of  poly¬ 
carbonate  and  extended  at  the  same  displacement  rate  of  8-467  X  10“  *  m/s. 
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Fig,  2.  Cooling/heating  of  polycarbonate  in  uniaxial  extension  at  a  displacement  rate  of 

8-467  X  10”’  m/s. 


the  resulting  temperature  as  a  function  of  time  would  be  very  different.** 
The  results  are  displayed  in  Fig.  2  in  which  cooling  terminates  at 
approximately  12  s  and  the  heating  portion  of  the  temperature  curve 
attained  only  a  gradual  rise  before  increasing  rapidly.  There  is  a  delay  in 
energy  dissipation  that  is  in  contrast  to  the  results  in  Fig.  1  for  the 
aluminum.  The  time  interval  associated  with  cooling/heating  is  sensitive  to 
changes  in  the  molecular  structure.  This  is  illustrated  by  the  difference  of  the 
two  curves  labelled  specimen  A  and  B  in  Fig.  2.  The  temperature  fluctuation 
data  also  monitor  the  change  in  the  chemical  composition  of  the  hardener 
molecules  and  the  curing  conditions  such  that  the  network  can  be  altered  in 
predictable  ways,  both  experimentally  and  theoretically.  Such  a  capability 
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is  made  possible  only  by  the  energy  density  theory®  and  can  be  used  to 
optimize  the  selection  of  the  fiber  and  matrix  material. 

2.3.  Interface 

Because  processing  can  vary  the  conditions  at  the  interface,  analytical 
modelling  would  not  be  complete  without  addressing  both  the  electro¬ 
chemical  and  mechanical  bonding  properties.  (The  assumption  that 
displacements  and  stresses  are  continuous  across  the  fiber  and  matrix  is 
obviously  unrealistic  and  leads  to  errors  that  depend  on  the  type  of 
loading. The  former  arises  by  wetting  the  solid  surface  with  a  fluid. 
Bonding  of  the  two  phases  is  achieved  through  intermolecular  forces. 
Theoretical  strengths  of  the  order  of  10®-10^MPa  can  result  within 
interatomic  distances  of  a  few  Angstroms  as  shown  schematically  in  Fig.  3. 
(Actual  strength  of  an  interfacial  bond  is  lower  because  the  classical 
theoretical  treatment  is  overly  idealized  and  is  not  able  to  include  the 
combined  effects  of  material,  geometry  and  load  transfer.)  The  Van  der 
Waals  bonds  are  slightly  weaker  with  theoretical  strengths  of  the  order  of 
10^-10^  MPa  and  involve  interatomic  distances  of  3-5  A.  Mechanical 
bonds  are  developed  as  a  result  of  differential  thermal  expansion  of  the  fiber 
and  matrix.  A  compressive  environment  is  created  that  results  in  high 
frictional  adhesion.  Compressive  loading  on  an  isolated  filament  in  a  glass 
fiber/epoxy  composite  can  result  in  a  frictional  bond  of  10®  Pa.  This, 
however,  is  spread  over  a  much  larger  distance  that  is  4-5  orders  of 
magnitude  larger  than  the  chemical  bond  (Fig.  3).  Minute  defects  such  as 
microvoids,  microcracks,  debonding,  etc.,  can  also  prevail  at  the  interface 
that  further  complicate  the  state  of  affairs  at  the  fiber/matrix  interface.  Even 
though  the  initial  interfacial  energy  states  arising  from  electrochemical  and 


Fig.  3.  Schematic  of  electrochemical  and  mechanical  interface. 
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mechanical  bonding  are  essential,  they  alone  are  not  sufficient  for 
characterizing  the  interface  integrity  and/or  failure.  The  true  nature  of  the 
interphase  depends  simultaneously  on  the  mechanics  at  a  two-phase 
boundary  and  rates  of  energy  transfer  that  are  application-specific.  No 
generalization  is  possible.  Each  case  must  be  analysed  separately  according 
to  the  specificity  of  the  system. 

The  consideration  of  load  transfer  across  interfaces  becomes  a  major 
issue  for  composites  in  structural  applications.  Unless  the  relations  between 
the  overall  bulk  behavior  of  the  composites  and  alleged  interfacial  response 
are  known,  no  design  guidelines  can  be  laid  down.  Depending  on  the  load 
type,’’  interface  properties  will  affect  the  composite  behavior  in  different 
ways.  This  is  illustrated  in  Fig.  4(a)  and  (b)  in  which  the  interfacial  energy 


Compression 


(a)  Compressive  lood 


(b)  Shear  load 


Fig.  4.  Interfacial  oscillation  in  energy  due  to  loading:  compression  and  shear. 
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oscillated  much  more  widely  in  shear  than  that  in  compression.  The  local 
material  elements  respond  at  rates  dictated  by  the  properties  of  the  two 
adjoining  phases  and  the  morphology  in  the  region  near  the  interface. 
Hence,  the  mutual  synergistic  effects  must  be  regarded  as  unknowns  and 
determined  theoretically.  The  fundamental  understanding  on  this  subject 
has  not  been  well  developed.  Apparent  contradictions  and/or  inaccuracies 
can  result  if  the  interface  properties  are  assumed  to  be  known  as  a  priori. 
This  occurs  even  in  the  case  of  a  single  phase  material  where  the  same 
constitutive  relation  is  assumed  to  prevail  at  each  point  of  the  system. 
(Conventional  mechanics  theories  assume  not  only  that  constitutive 
relations  are  known  but  they  are  the  same  at  each  point  in  the  system.  This 
invalidates  their  application  to  problems  where  the  local  strain  rates 
undergo  large  changes.)  The  stress  and  strain  relations  at  elements  near 
a  boundary  can  differ  widely  from  those  inside  a  body.®  They  must  be 
determined  separately  for  each  element  and  each  time  step  of  loading.  Such 
an  approach  must  be  adopted  for  analyzing  composite  interface  behavior. 

Debonding  of  fiber/matrix  interface  should  be  distinguished  from  that  of 
delamination  between  the  laminae  in  a  layered  composite  structure  (Fig.  5). 
The  amplitude  of  energy  oscillation  in  delamination  is  considerably  larger 
than  that  associated  with  debonding.  Since  both  normal  and  shear  stress 
may  both  be  present,  delamination  can  be  best  analysed  by  application  of 
the  energy  density  function  dlF/dK  The  stationary  values  of  dWjdV 
consider  failure  by  both  excessive  change  in  shape  and/or  fracture  as  they 
simultaneously  account  for  dilatation  and  distortion. 
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3.  DAMAGE  OF  COMPOSITES 

Composite  failure  can  be  best  addressed  in  terms  of  damage  over  an  area  or 
volume.  Above  all,  the  properties  of  the  constituents  as  measured  separately 
must  be  correlated  with  the  structural  features  of  the  composite.  This 
requires  a  precise  knowledge  of  how  the  time  dependent  energy  is 
distributed  and  used  to  deform  and  damage  the  composite  system. 

3.1.  Failure  Mechanisms 

Unlike  the  metals  where  failure  may  be  dominated  by  the  growth  of 
macrocracks,  fiber  reinforced  composites  fail  in  a  cumulative  fashion  that 
may  involve  a  combination  of  different  modes  such  as  fiber  breaking, 
matrix  cracking  and  fiber/matrix  interface  debonding  (Fig.  6).  Each  of  the 


Oebonding  Fiber 

breaking 


Matrix 

cracking 


Fig.  6.  Failure  modes  in  unidirectional  fiber  reinforced  composite. 


modes  may  occur  at  a  different  time.  Hence,  the  energy  associated  with 
creating  new  surfaces  in  a  unit  volume  at  a  given  instance  may  consist  of 


in  which  the  subscripts  f,  m  and  f/m  stand,  respectively,  for  the  fiber,  matrix 
and  fiber/matrix  interface.  (Other  forms  of  surface  energy  associated  with 
damage  may  be  involved  such  as  (dlF/dA),  that  arises  due  to  delamination 
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or  breaking  of  the  adhesive  in  a  layered  composite.)  Any  one  of  the 
component  surface  energy  densities,  say  {AW ! A may  be  zero  within 
a  certain  time  interval  if  debonding  does  not  occur.  In  other  words,  the 
proportion  of  fiber  breaking,  matrix  cracking  and  fiber/matrix  debonding 
changes  with  time.  What  is  available  to  create  new  damage  surface  at  the 
next  time  increment  is  {AW /A  A)*.  Hence,  the  total  surface  energy  density  is 


The  damage  accumulated  in  a  unit  volume  of  composite  element  can  be 
obtained  from 


(3) 


with  dIk'/dV'  being  the  area  under  the  true  stress  and  true  strain  curve.  The 
rate  change  of  volume  with  surface,  {AV/AA)^,  is  made  proportional  to  the 
respective  slopes  of  the  stress  and  strain  curves.  In  this  way,  the  energy  state 
for  each  composite  element  can  be  uniquely  determined  and  identified  with 
data  that  are  obtainable  from  uniaxial  tests.  This  is  accomplished  by 
referring  the  components  (d  ff'/dA)^,  (d  W/AA)^  and  (d  W/AA).  in  Fig.  7  to  the 


Fig.  7.  Schematic  of  a  unit  volume  of  composite  element. 
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damage  planes  such  that 


The  parameters  A,  B  and  C  are  determined  by  performing  three  different 
uniaxial  tests  along  the  axes  of  symmetry  for  the  composite  system.  In  this 
way,  the  constitutive  relation  for  each  composite  element  will  be  derived 
rather  than  assumed  as  it  is  now  being  done  by  the  classical  approach. 

Thresholds  for  characterizing  the  integrity  of  composite  systems  can  thus 
be  established  by  considering  the  rate  at  which  the  volume  and  surface 
energy  density  approach  their  respective  critical  values.  Because  of  eqn.  (4), 
it  suffices  to  consider  (d  WIAA)^.  Failure  of  composite  elements  is  assumed 
when 


and/or 


Depending  on  the  rate  of  load  transfer  in  the  composite,  the  exchange  of 
surface  and  volume  energy  density  can  vary  over  a  wide  range  as  indicated 
in  eqns  (5)  and  (6)  so  that  failure  mechanisms  involving  fiber  breaking, 
matrix  cracking,  interface  debonding,  etc.,  can  be  quantified  in  terms  of 
temperature  change,  energy  dissipation  and  the  rate  change  of  volume  with 
surface. 

Related  to  {dWIdA)^  and  (dIF/d/4)*  are  the  quantities  (dWIdV)^  and 
(dW/dV)*  which  are,  respectively,  the  dissipated  and  available  volume 
energy  as  defined  in  Fig.  8  such  that 


Note  that  (d W/d V)^  or  ^  is  the  area  oypq  and  (dW Id V)*  or  .s/  the  area  pgq. 
The  path  of  unloading  is  determined  analytically  in  the  energy  density 
theory.  It  is  of  interest  to  point  out  that  energy  dissipation  in  composite 
specimens  has  been  measured  experimentally.**  (One  of  the  experimental 
curves  for  ^  in  Ref.  16  became  negative  which  is  a  physical  impossibility 
and  obviously  caused  by  error  in  measurement.)  It,  of  course,  must  be 
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Fig.  8.  Schematic  of  true  stress  and  true  strain. 


positive  definite  and  increase  monotonically  with  time,  i.e. 

9^0;  — >0  (8) 

dt 

The  method  of  Electromagnetic  Discharge  Imaging  (EDI)  can  also  be  used. 
Once  9  is  known,  the  temperature  0  can  be  obtained  from 

A0  ,  AK  As 

— -  (9) 

0  AA  A9IAt 

in  which  Ae  is  the  increment  of  strain.  It  can  be  referred  to  the  plane  of 
homogeneity  (this  is  the  plane  defined  in  Fig.  7  such  that  (dVF/d4),  for 
i  =  ^,ri,C  are  related  as  shown  in  eqn.  (4))  with  Ae^  being  the  equivalent 
uniaxial  incremental  strain  for  a  multiaxial  stress  or  strain  state. 

Equation  (9)  has  been  used  successfully  for  determining  many  of  the 
previously  unexplained  thermal/mechanical  interaction  effects^' in 
isotropic  and  homogeneous  materials.  It  applies  equally  well  to  composites. 
The  quantity  0  is  not  equal  to  the  temperature  T  defined  by  classical 
thermodynamics.  It  is  equal  to  T only  when  all  energies  are  dissipated  in  the 
form  of  heat  Q.  In  general,  9  is  not  equal  to  Q.  The  strain  rate  of  energy 
dissipation  A9/Ae  can  be  associated  with  the  latent  heat  in  classical 
thermodynamics  for  determining  phase  transformation.  Again,  only  when 
9-*Q  will  the  following  relation  hold; 
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Here,  AQ  is  the  heat  required  per  unit  mass  of  substance  in  changing  the 
thermodynamic  state  from  the  temperature  T  to  T+ AT  and  Av  is  the 
change  of  specific  volume.  The  corresponding  change  in  entropy  is  Ay. 
Phase  transformation  is  a  rate  dependent  process,  an  effect  not  considered 
in  classical  thermodynamics. 

3.2.  Characterization 

A  unique  characterization  of  microfailure  mechanisms  in  composites 
can  be  made  by  specifying  the  quantities  AV/AA,  A©  and  A3,  all  of  which 
can  be  measured  individually  and  independently.  Their  combination  can 
uniquely  determine  any  changes  in  the  composite  microstructure.  There 
are  no  difficulties  at  present  in  measuring  AV/AA  via  a  knowledge  of 
the  displacement  field  and  local  temperature  changes  within  ±10"^  to 
+  10“^K  in  microelements  having  linear  dimensions  of  10'^  —  10“^  cm. 
The  response  time  of  the  electronic  instruments  must  be  adjusted  in 
accordance  with  the  loading  rate  that  governs  the  way  with  which  the 
material  microstructure  reacts.  A  major  challenge  is  the  detection  of  A3. 
Being  developed  at  the  Lehigh  Institute  of  Fracture  and  Solid  Mechanics  is 
the  electromagnetic  discharge  imaging  (EDI)  technique.*^  It  involves 
accelerating  the  electrons  surrounding  the  object  or  specimen  by  an 
external  field  such  that  the  air  molecules  are  ionized  to  enhance  an 
exponential  growth  in  the  number  of  electrons  and  positive  ions  causing  an 
avalanche. 

To  fix  ideas,  consider  the  composite  system  in  Fig.  9  consisting  of  three 
different  types  of  microfailure,  namely,  fiber  breaking,  matrix  cracking  and 


Fig.  9.  Hypothetical  microfailure  in  fiber  reinforced  composite. 
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fiber/matrix  debonding  in  addition  to  the  texture  of  fiber  and  matrix. 
A  sufficient  number  of  fibers  is  assumed  to  contain  in  the  macroelement  of, 
say  Ax  =  Ay  10“  ‘  cm,  such  that  the  microfailure  response  can  be  reflected 
by  the  time-dependent  macroscopic  parameters.  The  state  of  affairs  at  the 
microscopic  level  is  much  more  complex  because  microheterogeneities  are 
now  of  the  same  size  order  as  the  microelements,  say  Ax  =  Ay^  10“^  — 
10"^  cm.  Table  1  outlines  the  various  combinations.  The  variations  of 


Table  1 

A  hypothetical  characterization  of  microelements  for  a  fiber  reinforced  composite 

at  time  t^ 


Microelement  type 

(AF/AA)^ 

A© 

A^ 

Matrix 

Slightly  above 

Slightly  below 

Slightly  below 

Fiber 

Above 

Below 

Above 

Interface 

Slightly  below 

Slightly  above 

Near  average 

Fiber/matrix 

Near  average 

Near  average 

Slightly  above 

AF/AA,  A0  and  A^  for  the  four  different  types  of  microelements  can  be 
compared  with  the  average  values  obtained  by  considering  all  the 
microelements  in  the  macroelement.  The  loading  at  the  instance  tp  may  be 
such  that  AK/AA  and  Ai/  for  the  undamaged  composite  element  are  only 
slightly  above  the  average  while  A0  is  slightly  below.  This  is  illustrated  in 
Figs.  I0(a)-{c).  For  the  element  damaged  by  fiber  breaking  and  matrix 
cracking,  the  deviations  of  AF/AA,  A0  and  hsCC  from  the  average  may 
retain  the  same  trend  but  become  more  pronounced  as  shown  in  Figs.  1 1(a) 
-(c).  In  the  case  of  interface  debonding  in  Figs.  12(a)-(c),  AF/AA  maybe 
lower  and  A0  higher  than  the  average.  This  can  be  accompanied  by 


(a)  Volume  to  surface  chonqe  ( b) Temperoture  (c)  Dissipation 


Fig.  10.  Time  response  of  undamaged  composite  microelement. 
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0  tp  f  0  tp  t  0  tp  t 

(a)  Volume  to  surface  change  (b)  Temperature  (c)  Dissipation 


Fig.  11.  Time  response  of  composite  microelement  damaged  by  fiber  breaking  and  matrix 

cracking. 


0  tp  t  0  Ip  t  ®  *p  t 

(a)  Volume  to  surface  change  (b)  Temperature  (c)  Dissipation 


Fig.  12.  Time  response  of  composite  microelement  with  interface  debonding. 

a  higher  than  average  dissipation.  Each  microelement  will  exhibit  a  unique¬ 
ly  different  microstress  and  microstrain  response  even  for  the  same 
microstructure  detail  because  the  load  distribution  is,  in  general,  not 
uniform. 

Microstructure  change  can  also  take  place  by  phase  transformation.  This 
rate  process  depends  on  the  strain  rate  dissipation  energy  density  A^/Ae  in 
eqn.  (10)  at  a  given  tempierature  0.  Initiation  of  phase  transformation  can 
be  identified  with  the  reversal  of  curvature  of  the  JP-curve.  The  corres¬ 
ponding  time  f,  is  shown  in  Fig.  13(a).  The  values  of  (A^/Ae),  and  0^ 
at  f,  can  thus  be  obtained  from  the  curves  in  Fig.  1 3(b)  and  (c),  respectively. 
In  this  way,  the  distribution  and  change  of  the  different  phases  at  a  given 
location  can  be  assessed  analytically  and  experimentally. 

3J.  Scale  of  Observation 

The  transient  character  of  size/time/temperature  interaction  has  eluded 
those  working  in  composite  materials.  Failure  modes  observed  at  one  scale 


Advanced  Composite  Material  Behavior 


17 


(b)  Strain  rate  dissipation  (c)  Temperature 
energy  density 


Fig.  13.  Assessment  of  local  phase  transformation  from  the  Jt^-curve. 


level  may  differ  from  that  seen  at  another  level  and  they  are  further  com¬ 
plicated  by  changes  in  loading  rates.  In  general,  dilatation/distortion  and 
cooling/heating  tend  to  flip-flop  as  the  scale  level  of  observation  is  altered; 
depending,  of  course,  also  on  the  time  response.  The  fundamentals  of  this 
alternating  mechanism  are  discussed  in  Ref.  18  and  will  only  be  mentioned 
briefly  in  relation  to  an  element  dominated  by  macrodilatation  and  another 
by  macrodistortion  as  shown,  respectively,  by  Figs.  14(a)  and  (b).  (The 
combination  of  size/time/temperature  data  are  selected  arbitrarily  for  this 
discussion.  Actual  values  for  uniaxial  tensile  and  compressive  metal 
specimens  have  been  obtained  and  can  be  found  in  Ref.  19.)  These  two 
elements  are  necessarily  located  at  different  locations  of  the  same  system. 
To  be  emphasized  in  Fig.  14(a)  is  that  macrocooling  and  macrodilatation 
for  response  times  of  1-10  s  correspond  to  microheating  and  micro¬ 
distortion  when  the  same  event  is  viewed  within  the  time  interval  of 
10“^-10~‘s.  (Macro  refers  to  dimensions  of  10~^-10"^cm;  micro  to 
10“^-10"*cm;  and  atomic  to  10~*-10”®cm.  The  size/time/temperature 
scale  can  shift  depending  on  the  loading  or  local  strain  rate.)  Disturbances 
at  the  atomic  scale  refer  to  even  smaller  response  time  and  temperature 
fluctuation.  The  situation  for  the  element  in  Fig.  14(b)  is  opposite  to  that  in 
Fig.  14(a).  What  was  macrocooling  now  becomes  macroheating.  The  same 
applies  to  the  micro-  and  atomic-scale  together  with  the  corresponding 
time  response  and  temperature.  Scaling  of  size/time/temperature  must  be 
assessed  quantitatively  such  that  seemingly  different  behavior  of  the  same 
failure  process  when  viewed  at  the  atomic,  microscopic  and  macroscopic 
level  can  be  related  uniquely  by  the  three  variables  mentioned  earlier. 
Microstructure  change  must  be  addressed  in  its  entirety  with  reference  to 
the  specimen  response.  Whether  its  influence  could  be  adequately  reflected 
by  the  macroscopic  variables  or  not  depends  on  the  specific  application. 
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Atomic 

cooling- 

Atomic 

dilototion 

(iO"*°K) 


.  Ambient 
Condition 
Mocrccooling:  Mocrodilototion 
(lO'^'K) 


Response  time~  1  to  10  sec. 
Response  time= 

Response  time=  10  '‘sec. or  smaller 


(o)  Response  of  element  olong  path  of  crock  growth 


(b)  Response  of  element  olong  direction  of  moximum  shope  change 

Fig.  14.  Size/time/temperature  response  of  element  at  different  scale  level,  (a)  dominated  by 
macrodilatation  and  (b)  dominated  by  macrodistortion. 


4.  ELECTROMAGNETIC  DISCHARGE  IMAGING: 
NONDESTRUCTIVE  EVALUATION 

Because  of  the  complex  nature  of  composite  failure  behavior,  a  variety  of 
analytical  models  have  emerged  in  an  attempt  to  explain  the  different 
failure  modes.  Since  no  unified  treatment  has  come  forth,  damage 
inspection  procedures  have  relied  solely  on  the  development  of  hardwares 
with  the  assumption  that  the  data  related  to  debonding,  fiber  and  matrix 
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degradation,  etc.,  will  eventually  be  understood  and  analysed  by  theories. 
Among  the  techniques  used  are  ultrasonics,  radiography,  acoustic  emis¬ 
sion,  thermogfaphy,  eddy  current  and  a  host  of  others.  Such  a  procedure 
has  been  loosely  referred  to  as  Non-Destructive  Testing  (NDT).  What  is 
urgently  needed  is,  of  course,  a  consistent  way  of  assessing  composite 
damage  so  that  the  seemingly  different  failure  modes  can  be  weighed  on 
a  common  basis.  There  is  no  other  way  for  qualifying  the  integrity  of 
composite  materials  and  structures. 

Damage  being  a  time-rate  dependent  process  is  more  widely  spread 
throughout  the  composite  simply  because  the  load  reacts  more  sensitively 
with  material  inhomogeneity.  The  idea  of  NDT  is  to  correlate  the  difference 
between  the  emitted  and  received  signal  with  the  size  and  location  of 
damage.  Those  energies  that  are  already  present  within  the  system  may  also 
be  released  during  the  damage  process  to  interrupt  the  NDT  signal.  The 
unique  identification  of  energy  dissipation  rate  with  failure  mode  is 
a  prerequisite  for  the  quantitative  assessment  of  damage,  i.e.  the  quantity 
dS/dr  or  those  in  eqn.  (1)  involving  dfdW'/d/llf/dt,  d{dW ldA)Jdt,'etc.  The 
Electromagnetic  Discharge  Image  (EDI)  technique commonly  known 
as  Kirlian  photography,  is  particularly  suited  for  incorporating  the  energy 
density  theory.®  (This  technique  has  been  used  extensively  in  the  USSR  in 
NDT  and  is  rarely  publicized.  Preliminary  efforts  made  by  the  US  Navy 
and  Lawrence  Livermore  Laboratory  lacked  the  theoretical  support  in 
NDT  application.)  Utilized  in  the  process  is  a  high-voltage  electrical 
discharge  field.  Electrons  with  the  specimen  are  accelerated  and  multiply 
exponentially.  Streamers  are  thus  created  that  move  at  speeds  of  10’  — 10® 
cm/s.  In  air,  at  high  field  strength,  the  streamers  are  blue  while  a  reddish- 
purple  glow  appears  in  low  eledtric  fields.  The  intensity,  color  and  pattern  of 
the  discharge  image  can  be  recorded  photographically  or  digitally  and 
identified  with  the  mechanical,  thermal,  chemical  and  electrical  disturb¬ 
ances  in  composites  as  they  are  damaged.  Mechanical  discontinuities  in  the 
range  of  10~  ‘-10~®  mm  can  be  detected.  Such  a  wide  range  of  detection 
capability  is  unmatched  by  the  other  conventional  methods.  The  flow  chart 
in  Fig.  1 S  shows  how  the  available  (or  dissipated)  energy  is  monitored  by 
EDI.  A  laboratory  EDI  model  has  been  designed  by  the  Institute  of 
Fracture  and  Solid  Mechanics  at  Lehigh  University.  The  combined 
thermal,  mechanical,  chemical  and  electrical  effects  arc  recorded  every¬ 
where  on  a  specimen  as  indicated  in  Fig.  16(a).  Spectral  analysis  of  the  pixel 
can  then  be  made  to  obtain  the  optical  intensity  as  a  function  of  radiation 
wavelength  (Fig.  16(b)).  For  each  radiation  wavelength  the  optical 
intensity  /j  can  be  isolated  and  identified  with  a  given  energy  level  Uj.  This 


Electrical 


Cause  space  

Fig.  15.  Cause  and  effect  on  intensity  and  color  of  emitted  radiation. 


(a)  Combined  intensity  (b)  Individual  wavelength 

Fig.  16.  Intensity  versus  radiation  wavelength. 
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is  accomplished  by  means  of  a  plasma  spectral  analysis,  the  results  of  which 
are  shown  schematically  in  Fig.  17(a)-(c).  The  results  can  be  further 
processed  and  displayed  by  contour  plots  in  two  or  three  dimensions  as 
illustrated  in  Fig.  18  by  which  the  individual  thermal,  mechanical,  chemical 


(a)  Energy  (b)  Energy  (c)  Energy 


Fig.  17.  Variations  of  intensity  with  energy  level. 


or  electrical  effect  can  be  sorted  out.  The  color,  intensity  and  distribution  of 
the  contours  provide  information  on  the  location  and  magnitude  of  the 
stationary  values  of  the  energy  density  whose  relation  with  material 
damage  has  been  discussed  earlier.  These  contours  will  change  as  a  function 
of  time,  if  load  is  present  and  hence  the  influence  of  time-rate  dependency  on 
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damage  can  also  be  evaluated.  To  be  more  specific,  oscillation  of  the 
electrochemical  bond  energy  may  be  displayed  by  contour  plots  in 
a  two-dimensional  space.  Refer  to  Fig.  19  in  which  contours  with  high  and 
low  energy  densities  represent  fluctuations  in  the  electrochemical  forces. 
They  can  be  measured  by  the  EDI  technique  and  calculated  by  the  energy 
density  theory.^  As  the  scale  of  investigation  is  enlarged  near  the 
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Fig.  19.  Schematic  of  electrochemical  interface  energy  density  contours  detected  by  EDI. 


fiber/matrix  interface,  stationary  values  of  the  mechanical  and  thermal 
energy  density  can  be  obtained  in  the  same  way.  The  same  theoretical  and 
experimental  method  must  be  applicable  to  cover  the  range  of  energy 
dissipation  rates  that  occur  in  composite  material  damage. 
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ABSTRACT 

An  analysis  of  failure  of  composites  originating  from  rigid  inclusions  is 
undertaken.  The  following  cases  are  considered:  ( i )  inclusions  with  cuspidal 
points  embedded  in  an  elastic  matrix  and  ( ii )  inclusions  partially  bonded  to 
the  matrix.  In  the  first  case  failure  initiates  from  the  cuspidal  points,  while  in 
the  second  case  from  the  tips  of  the  interfacial  crack  which  coincides  with  the 
unbonded  part  of  the  inclusion  due  to  the  existing  stress  singularity  at  these 
points.  A  closed  form  solution  for  the  problem  of  partially  bonded  inclusions 
is  obtained  by  using  the  conformal  mapping  technique  of  the  complex  variable 
theory  of  elasticity  and  reducing  it  to  a  Hilbert  problem.  Results  are  obtained 
for  special  inclusion  shapes  including  the  fiber,  the  hypocycloidal,  the 
astroidal,  the  square  and  the  triangular  inclusion.  The  strain  energy  density 
theory  is  used  to  study  the  failure  of  the  composite  and  the  critical  load  and 
initial  fracture  angles  are  determined. 


1.  INTRODUCTION 

Composites  made  up  of  a  soft  matrix  and  stiff  particles  are  analysed 
regarding  their  mechanical  behavior  along  two  distinct  approaches:  the 
macro-approach  and  the  micro-approach.  The  macro-approach  considers 
the  composite  as  a  homogeneous  system  which  has  the  overall  material 
properties  of  the  macrostructure  and  applies  the  principles  of  continuum 
mechanics  of  homogeneous  (isotropic  or  anisotropic)  materials.  The 
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micro-approach  takes  into  account  the  detailed  geometrical  and  physical 
properties  of  the  constituent  materials  and  models  the  real  system  by 
a  simpler  one.  The  approach  chosen  in  the  analysis  of  composites  depends 
on  the  nature  of  the  problem  and  the  objective  of  the  study. 

In  the  analysis  of  composites  within  the  framework  of  the  micro-approach 
certain  idealizations  referred  to  the  geometry  and  the  properties  of  the 
constituent  materials  are  usually  made.  Thus,  the  filler  particles  are 
simulated  by  bodies  with  simple  geometrical  shape,  like  spheroidal, 
ellipsoidal  or  cylindrical.  Furthermore,  when  their  modulus  of  elasticity  is 
much  higher  than  the  modulus  of  elasticity  of  the  matrix  the  particles  are 
considered  stiff. 

In  many  composites  the  reinforcing  constituents  are  of  irregular  shape 
with  sharp  angles,  like  the  various  inorganic  fillers,  the  metal  or  boron 
filaments,  the  aggregate  or  sand  particles  in  concrete.  In  such  cases  in  the 
sharp  angle  corners,  high  stress  concentrations  develop  and  therefore  they 
are  nuclei  for  the  generation  of  cracks  and  slip  bands  leading  to  failure. 
A  commonly  observed  failure  mode  of  multiphase  uiaterials  is  the 
debonding  of  the  different  phases  due  to  manufacturing  and/or  loading 
conditions,  thus  forming  interfacial  cracks.  This  is  usually  encountered  in 
concrete  where  cracks  are  formed  along  the  boundaries  of  the  aggregates 
which  are  embedded  in  the  mortar. 

In  the  present  communication  in  an  attempt  to  model  the  failure 
behavior  of  composites  originating  from  rigid  inclusions  within  the 
framework  of  micro-approach  the  following  two  problems  arc  studied; 

(i)  Rigid  inclusions  with  cuspidal  points  along  their  boundaries  embed¬ 
ded  in  an  elastic  matrix.  In  such  cases  high  stress  concentrations  are 
developed  in  the  vicinity  of  the  cuspidal  points  which  constitute 
nuclei  for  failure  initiation. 

(ii)  Rigid  inclusions  of  general  shape  partially  bonded  to  an  elastic 
matrix.  The  unbonded  part  of  the  inclusion  constitutes  an  interfacial 
crack  from  which  failure  of  the  composite  starts. 

In  both  cases  the  stress  field  in  the  composite  is  first  analysed  and  it  is 
then  coupled  with  an  appropriate  failure  criterion.  Due  to  the  importance 
of  the  problem  elasticity  solutions  of  a  number  of  geometrical  configur¬ 
ations  of  inclusions  perfectly  or  partially  bonded  to  a  matrix  have  appeared 
in  the  literature.  Panasyuk  et  n/.‘  studied  the  case  of  a  rigid  inclusion  with 
cuspidal  points  on  its  boundary  embedded  in  an  elastic  matrix.  They  found 
that  the  stress  field  in  the  vicinity  of  the  cuspidal  points  presents  an  inverse 
square  root  singularity  and  it  is  expressed  in  terms  of  two  stress 
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concentration  factors.  The  problem  of  a  rigid  inclusion  partially  bonded  to 
an  elastic  matrix  has  been  studied  by  various  researchers.^  Circular, 
elliptical  and  general  shaped  inclusions  have  been  considered  and  attention 
has  been  paid  to  the  neighborhood  of  the  ends  of  the  unbonded  part  of  the 
inclusion. 

In  this  work  the  failure  of  composite  plates  containing  a  perfectly  or 
partially  bonded  rigid  inclusion  in  an  elastic  matrix  is  analysed.  A  general 
solution  for  the  stress  and  displacement  field  is  obtained  for  a  curvilinear 
inclusion  with  an  interfacial  crack  using  the  method  of  conformal  mapping 
in  conjunction  with  the  analytic  continuation  of  the  complex  potentials. 
The  problem  is  reduced  to  a  Hilbert  problem  for  one  of  the  complex 
potentials  and  general  formulae  for  the  determination  of  the  various 
unknown  coefficients  of  the  solution  are  given.  Special  cases  include  the 
square  and  the  triangular  inclusions.  The  results  of  the  stress  analysis  are 
combined  with  the  strain  energy  density  failure  criterion.**  *^  The  critical 
load  for  fracture  initiation  and  the  more  vulnerable  failure  sites  are 
determined.  The  failure  characteristics  of  the  composite  as  influenced  by  the 
type  and  geometrical  configuration  of  the  inclusion,  the  material  properties 
and  the  type  of  loading  are  disclosed. 


2.  RIGID  INCLUSIONS  WITH  CUSPIDAL  POINTS 


2.1.  The  Stress  Field 

Consider  a  rigid  inclusion  with  a  cuspidal  corner  O  perfectly  bonded  to 
an  infinite  isotropic  elastic  plate  which  is  subjected  to  a  system  of  stresses  at 
infinity  (Fig.  1).  A  reference  frame  of  Cartesian  coordinates  is  attached  to 
the  point  O  with  the  x-axis  along  the  tangent  of  the  boundary  of  the  in¬ 
clusion  at  O.  The  local  stress  field  in  the  vicinity  of  the  point  O  is  given  by;* 
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where  the  coefficients  k^  and  ^2  ate  independent  of  the  coordinates  r,  0  and 
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Fig.  1.  Geometry  of  a  general  shaped  inclusion  with  a  singular  corner  O  embedded  in 

a  matrix. 


depend  only  on  the  loading  conditions,  the  material  of  the  plate  and  the 
geometrical  shape  of  the  inclusion  at  the  cuspidal  corner.  In  the  above 
relations,  k  =  (3  — v)/(l  +  v)  or  k  =  3— 4v  for  plane  stress  or  plane  strain 
conditions  respectively,  with  v  representing  the  Poisson’s  ratio  of  the 
material  of  the  plate.  Note  the  inverse  square  root  singularity  of  the  stress 
field  near  the  point  O. 

Failure  of  the  composite  plate  initiates  from  the  cuspidal  point  O  due  to 
the  high  intensification  of  the  stress  field  around  this  point.  For  the 
determination  of  the  critical  failure  loads  and  the  fracture  path  use  is  made 
of  the  minimum  strain  energy  density  criterion. 

2.2.  The  Minimum  Strain  Energy  Density  Criterion 

Failure  initiation  from  the  cuspidal  point  of  the  inclusion  is  described  by 
the  minimum  strain  energy  density  criterion  proposed  by  Sih  “  •  ‘  ^  and  used 
for  the  solution  of  a  host  of  mixed-mode  crack  problems  by  Gdoutos.  *  ^ 

The  fundamental  quantity  for  unstable  fracture  is  the  strain  energy 
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density  factor  S  defined  by 


(2) 


where  dW^dF  is  the  strain  energy  density  per  unit  volume  and  r  is  the 
distance  from  the  cuspidal  point.  For  the  plane  elastic  problem  dH^dFis 
expressed  by 


dW 

dF 
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where  /x  represents  the  shear  modulus. 

It  is  assumed  that  the  fracture  path  from  the  cuspidal  point  of  the 
inclusion  follows  the  direction  of  minimum  strain  energy  density  factor, 
defined  by 


d^S 

dd^ 


>  0 


(4) 


Unstable  fracture  occurs  when  the  minimum  value  S„j„  of  S  becomes 
equal  to  a  critical  value  S„  which  is  a  material  constant,  that  is 

=  (5) 

Relations  (4)  define  the  fracture  path  originating  from  the  cuspidal  point, 
while  eq.  (5)  gives  the  critical  load. 

From  eqns  (l)-(3)  is  obtained 

S  =  aijfci +20,2^1^2  +  022^2  (6) 

where  the  coefficients  o.^- (1,7  =  1,2)  are  given  by: 

0 

Ib/xo,,  =2(k:  —  1)cos^-+k:^  +  (2)c+  1)cos^  6 

l6/xai2  =  —[(K—l)+2KCOS0]sin0  (7) 

0 

16/xa22  =  2(>c—  l)sin^-+K^— (2k—  l)cos^  0 


The  minimum  strain  energy  density  criterion  is  used  for  the  deter¬ 
mination  of  fracture  of  a  composite  plate  with  a  linear,  an  astroidal  and 
a  hypocycloidal  inclusion. 


2J.  The  Linear  Inclusion 

A  rigid  rectilinear  inclusion  of  length  2/  is  embedded  in  an  elastic  plate 
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which  is  subjected  to  a  uniform  uniaxial  stress  a  making  an  angle  ^  with 
the  axis  of  the  inclusion.  The  stress  concentration  factors  and  ^2  are  given 
by  ‘ 


P— 

Fig.  2.  Variation  of  the  normalized  critical  stress  of  fracture  for  the  case  of 

a  fiber  inclusion  versus  its  orientation  angle  p  according  to  the  minimum  strain  energy  density 

criterion. 
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From  eqns  (4)  -  (7)  the  fracture  path  and  the  critical  stress  ff„  for  unstable 
fracture  of  the  composite  plate  are  obtained.  Figure  2  presents  the  variation 
of  the  normalized  critical  stress  versus  the  inclination  angle  of  the 
inclusion  P  for  various  values  of  the  material  constant  k.  Note  that  for  jc  =  1, 
(T„  is  symmetrical  with  respect  to  P= 45°  for  which  becomes  infinite.  For 
1C  =  1-4  and  1-8,  becomes  maximum  in  the  interval  0<^<90,  while  for 
the  remaining  values  of  k  it  takes  its  maximum  value  for  P  =  90°.  The 
variation  of  the  initial  fracture  angle  is  shown  in  Fig.  3. 


Fig.  3.  Variation  of  the  fracture  angle  0„  versus  the  angle  P  of  the  orientation  of  the  fiber 
according  to  the  minimum  strain  energy  density  criterion. 
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2.4.  The  Astroidal  Inclusion 

An  astroidal  inclusion  is  circumscribed  by  a  circle  of  radius  a  and  it  is 
referred  to  a  reference  frame  Oxy  with  O  coinciding  with  the  center  of  the 
circle  and  the  singular  comer  j=0  lying  on  the  x-axis.  The  equation  of  the 
inclusion  with  respect  to  the  frame  Oxy  is  given  by 

with  2  =  re'*  and  (=e'*. 

If  P  is  the  angle  that  the  applied  stress  tr  subtends  with  the  x-axis,  then  the 
ki  and  ^2  stress  concentration  factors  are  given  by‘ 


1  3*:  ,  . 


*2 ’■  ^12  j^sin  W-2^) 


with  j  =  0,  1,  2,  3  for  the  four  cuspidal  points  of  the  inclusion. 

Because  of  the  existing  symmetry  only  the  critical  fracture  stresses 
from  the  points  j=0  and  j  =  1  are  determined.  Fracture  of  the  composite 
plate  starts  from  the  more  vulnerable  cuspidal  point  which  has  the  smaller 
critical  stress.  Figure  4  presents  the  variation  of  the  dimensionless  quantity 
ff„(3a/256/iS„)‘'^  versus  angle  p.  Regions  where  fracture  starts  from  either 
of  the  cuspidal  points  y = 0  or  7  =  1  are  separated  by  a  dotted  line. 

2.5.  The  Hypocycloidal  Inclusion 

The  hypocycloidal  inclusion  circumscribed  by  a  circle  of  radius  a  and 
referred  to  a  reference  frame  Oxy  with  O  coinciding  with  the  center  of  the 
circle  and  the  Ox  axis  being  the  axis  of  symmetry  of  the  inclusion  is 
described  by  the  following  equation. 


with  2  =  re'*  and  C  =  e‘®. 

If  P  is  the  angle  that  the  applied  stress  a  subtends  with  the  x-axis,  then  the 
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Fig.  4.  Variation  of  the  normalized  critical  stress  of  fracture  o„(3a/256ftSc,)‘'^  versus  angle 
P  for  the  case  of  an  astroidal  inclusion  according  to  the  minimum  strain  energy  density 
criterion.  Regions  where  fracture  starts  from  either  of  the  comers  j=0  or  j=  \  are  indicated. 


ki  and  fcj  stress  concentration  factors  are  given  by  ' 

with  j  =  0,  1,  2  for  the  three  cuspidal  points  of  the  inclusion. 
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As  in  the  previous  cases  the  critical  stress  and  the  fracture  angle  Og 
from  the  more  vulnerable  cuspidal  point  of  the  inclusion  are  determined. 
The  variation  of  the  quantity  versus  angle  ^  for  fracture 

initiation  from  the  cuspidal  points  j=0  and  ;  =  2  is  shown  in  Fig.  5. 


Fig.  5.  Variation  of  the  normalized  critical  stress  of  fracture  ir„(a/72fiS^,y  ^  versus  the 
orientation  angle  of  the  inclusion  ^  for  the  case  of  a  hypocycloidal  inclusion.  Regions  where 
fracture  starts  from  either  of  the  comers  7=0  or  j  =  2  are  indicated. 
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3.  PARTIALLY  BONDED  RIGID  INCLUSIONS 


3.1.  Statement  of  the  Problem  and  Hilbert  Formulation 

Consider  a  rigid  curvilinear  inclusion  perfectly  bonded  to  an  elastic 
infinite  matrix  except  from  a  part  As  of  its  boundary  forming  an  interfacial 
crack  (Fig.  6).  Denote  by  Aq  the  bonded  part  of  the  inclusion  boundary,  and 
by  K  and  fi  the  elastic  properties  of  the  matrix.  A  system  of  uniformly 
distributed  principal  stresses  and  is  applied  at  infinity  where  the 
direction  of  makes  an  angle  cp  with  the  x-axis. 


/ 


Fig.  6.  Geometry  of  a  rigid  curvilinear  inclusion  partially  bonded  to  an  elastic  matrix  and  its 
conformal  mapping  on  to  the  unit  circle. 

The  matrix  occupying  the  z-p!ane  is  mapped  in  the  infinite  region  Z  of  the 
C-plane  bounded  by  the  unit  circle  T  by  means  of  the  function  r  =  m(C) 
which  has  the  form 

z  =  R^C  +  y  +  ^+---+^^  (13) 

where  R  is  a  real  and  h,,  hj-  •  •  • ,  are  generally  complex  constants.  The 
values  of  these  constants  are  determined  so  that  the  contour  A  in  the 
2-plane  corresponds  to  the  circumference  T  of  the  unit  circle  in  the  C-plane. 
The  positive  sense  of  describing  the  contour  A  is  chosen  to  be  clockwise  so 
that  the  region  S  remains  on  the  left  when  moving  in  the  positive  direction. 
Putting 


^=e^(cos>?-l-isin  rj) 


(14) 
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circles  of  radii  and  straight  lines  tj  =  ri^  in  the  C-planc  define  an 

orthogonal  curvilinear  coordinate  system  (4,  rj)  in  the  z-plane.  Denote  by 
r,=rie‘*‘  and  t2  =  r2e‘®^  the  tips  of  the  interfacial  crack  A,  which  are 
mapped  on  the  unit  circle  F  at  points  ffi=e'*‘  and  <72=®'*^- 
Assuming  that  the  crack  lips  are  stress-free  the  boundary  conditions  of 
the  problem  take  the  form 

u(<T)-l-it)(ff)=iem(ff)  (reFu  (15) 

af^ar)+i<Tf^(ff)=0  (reFs  (16) 

where  u  and  v  are  the  Cartesian  components  of  the  displacement,  and 
denote  stress  components  refered  to  the  system  (ij,  tj)  and  e  represents  the 
rotation  of  the  inclusion. 

Using  the  complex-variable  formulation  of  the  plane  elastic  problem  for 
curvilinear  boundaries'^  eqns  (15)  and  (16)  give  the  following  equations 

lFo(ff)-lFort^)=0  ffeFs  (17) 

KWo{a)+W^a)=4ifiEm'(a)  asT^^  (18) 

for  the  complex  function  Wo(C)=m'(C)lF(0  where  is  the  usual  complex 
potential  of  the  theory  of  elasticity.  Wq  (a)  and  fF"  (<^)  denote  the  limit 
values  of  lFo(0  as  C  tends  to  <r  from  L  and  R  respectively,  where  L  represents 
the  image  region  of  the  matrix  and  R  the  circular  hole  in  the  C-plane. 
Equations  (17)  and  (18)  constitute  a  nonhomogeneous  Hilbert  problem 
with  line  of  discontinuity  the  arc  F^  described  in  a  clockwise  sense.  It  is 
obtained  for  the  unknown  function  (Ref.  14) 

Wo(0=^XmO  +  X{OR{Q  (19) 

fc 


where 


m'(C)d<T 

X‘-f<rX<T-C) 


R(C)=(A,C  +  A„)+^+^  +  -  +  ^^^ 


X(0= 


(C-ff2)^~‘ 


T=0-5-l-iA, 


A= 


lOgK! 

2n 


(20) 

(21) 

(22) 


(23) 
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The  function  X(0  is  the  Plemelj  function  of  the  problem  and  is 
holomorphic  in  the  whole  C-plane  cut  along  the  arc  To  on  which 
A'*((r)=  — /cA^‘'(<t).  Only  the  branch  of  the  function  XiQ  for  which 
limj^oo[CA'(0]  =  l  will  be  considered  in  the  sequel. 

In  calculating  the  integral  /(C)  it  is  observed  that  the  function 
F{Q=m’(0/X(0  is  holomorphic  in  the  whole  plane  cut  along  except  at 
the  points  C  =  oo  and  C=0  at  which  it  has  poles  of  orders  1  and  (p+ 1), 
respectively.  Thus,  the  principal  parts  of  F{0,  ffiiO  and  ^2(0  at  C  =  oo  and 
C=0  will  have  the  form 


g,(0=S,C  +  So 


92(0  =  ^  + 


(24) 

(25) 


Using  a  contour  A  surrounding  the  arc  TfllRef.  14)and  observing  that  on 
rDX*(<T)=  ~KX^(a)  the  following  is  obtained  for  the  integral  /(C) 


1(0= 


K  r  F(C)dff 

2jii(l+K:)J^  ff-C 


which  by  the  well-known  properties  of  the  Cauchy  integral  gives 


(26) 


m= 


1+K 


m'(C) 

2r(C) 


-gAO 


-92(0 


(27) 


Introducing  this  value  of  /(C)  into  eqn.  (19)  we  obtain  for  the  function 

W'o(C) 

Wo(0  =  <?[m'(C)  -  (g  i(C) + g2(0]2f  (C)]  +  K(C)X(C)  (28) 

where 


4i/ie 
1  -t-ic 


(29) 


Equation  (29)  gives  a  closed  form  solution  for  the  unknown  function 
Wq(Q  from  which  the  stress  and  displacement  field  are  determined. 


3.2.  Determination  of  the  Unknown  Coefficients 

The  unknown  coefficients  /I ,,  >1  _  i,  /4  _  .  ,  /4  _  ,  of  the  function 
R(C)  and  the  rotation  e  of  the  inclusion  are  determined  from  the  conditions 
that  the  complex  potentials  of  the  problem  should  be  holomorphic  in  the 
region  L  of  the  C-plane  and  have  a  particular  behavior  at  infinity.  After 
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lengthy  calculations  there  is  obtained 


^0  — 


2 


(30) 


where  is  the  rotation  at  infinity,  while  the  (p  + 1) coefficients  A.i,A_2, 
are  determined  from  the  following  (p+ 1)  equations 

5,-(p-l)/C,_i=0  (31) 


B,-2K2=  0 

B2-K,=R(N^-T^)c-^>^ 

5i=0 

iM-i-  N' 

Bo  +  K.,  = - - - Rb,{N^-T^)c-^^* 

where  M  is  the  moment  of  the  stresses  applied  on  A  about  the  origin,  N'  is 
a  real  constant,  and 


#>-m+  1 


a  =  0 

(32a) 

/  p+l  \  p+1 

Bo  =  MR-YS-A]+TA-sd, 

\  s=0  /  s=0 

(32b) 

p-m 

^m=  X  T^-»ap-„-,(m=-l,0,  1, ...,  p) 

s  =  0 

(33) 

'  II 

II 

1 

(34) 

7;-„  =  5',-„+  I^h*W.,,(m  =  2,  3,...,(p+l)) 

*=  I 


with 

h-i  =  l  bo=0 

r  m+l  T  m+l 


(35) 
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(0 

(1  — tU(Ti)'(l  — (Uffjf 

(p+l)Df(0)  DV”(0) 

"^‘■(p+l)!"  (p  +  1)!  “  p! 

D<^+i)(0)=  t 


(36a) 

(36b) 


>l(fl))  = 


xa 


1 


1— COffi  1—0)02 

X<*'(0) 

k\ 


A'(0)=e^<^*’“>"‘*'’ 


i^«')(0)  =  2(— l)*'^‘/c!r  cos 


ik+  I 


!g-i(*+l)«o 


r  =  yo-25  +  )S  =  tan  ‘  2A,  0^P<2ti 

w  =  0,-0i,  00  =  ^4^ 


(37) 

(38) 

(39a) 

(39b) 

(40) 

(41a) 

(41b) 


Furthermore,  the  coefficients  of  the  functions  p,(C)  and  P2(C)  defined  by 
eqns  (24)  and  (25)  are  given  by 


Si  =  R,  So=-RD2 


(42) 


and 


with 


lff=C 


c  =  lO-Rbi-2Rb2...  -pRb^y 


f=[S-,S_2 

...  s 

ir 

1 

O 

dt 

...  d,  - 

n  = 

0 

do 

...  d,., 

-0 

0 

do  - 

(43) 


(44) 
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3.3.  Local  Stress  Distribution 

The  analysis  of  the  stress  distribution  in  the  neighborhood  of  the  crack 
tip  is  of  particular  importance  and  allows  study  of  the  growth  character¬ 
istics  of  the  interfacial  crack.  Consider  the  tip  t2  of  the  crack  and  the  tangent 
t2X  of  the  inclusion  at  the  point  tj  (Fig.  7).  Using  the  asymptotic  expansion 


Fig.  7.  Crack  tip  coordinate  system. 


of  the  function  and  its  analytic  continuation  for  jCI  <  1  around  the 
point  tj  in  conjunction  with  the  equations  relating  %  (0  to  the  stress  field 
and  after  lengthy  algebra  the  following  equations  of  the  curvilinear  stress 
components  cr,,  and  are  obtained 

—  j^cos  ^  Ain  2n^ + e^^*"  "  cos  —  Ain  2np 
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I 


r . 


H - sin  -+Aln  2np 

2v^L  \2 


+ Ain  2npJ — *’  sin  — Ain  2np'j 


-  yj\  +4A^  sin  a  cos  + Ain  2np  — 

<T^^=— ^=1^3  cos  + Ain  27cp^— e^^‘’'““*cos^^— Ain  2np 

—  yj\  +4A^  sin  a  sin  +  Ain  2np  —  (p^  Ki 


(45) 


<T.  = 


-  — j^3  sin  +  Ain  2nf^ + e^^'"  “•  sin  —  Ain  2np 

-yj\+  4A^  sin  a  cos  +  Ain  2np  —  cp^  /Cj 

sin  + Ain  2kp^— e^'*‘’'''”sin  — Ain  27ip^ 


(46) 


A  Ad 


2y/^' 


+  ^1  +4A^  sin  a  cos  +  Ain  2np  -  (p 

e^“  r  -I  \ 

- P=  cos  -  +  /.In  2np  I 

2v^L  V2 

(3(x 

— +Aln  2np  —  (p 
<p  =  tan  "  '2A 

The  stress  intensity  factors  K,  and  are  given  by 


2 


+  e^^*’'‘'”cos(  -  -Ain  2np 


K, 


(47) 

(48) 


K ,  =  e - - — (A  COS9  +  B  sin 9) 
\lsmf 


(49) 


K,  =  e  - — (/!  sin  0  —  B  cos  0) 

^  -'sin^ 
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where 

0  =  0o+~^“+'^l*'j^47tsin^lm(ff2)lJ.  co-O^-Oi,  20o  =  0i  +  02 
Re[m'(<T2)]Re[P(ff2)]  +  lm[m'(g2)]  lm[P(g2)] 

Re[w'(ff2)]  Ini[P(g2)]-Im[m'(g2)]  Re[P(g2)] 

with 

P(C)  =  R(C)-^[S, (0  +  92(0] 

1  T  ?c 


(50) 


(51) 


3.4.  The  Square  Inclusion 

A  rigid  curvilinear  rounded-off  angle  square  inclusion  with  a  sym¬ 
metrically  located  interfacial  crack  is  embedded  in  a  plate  subjected  to 
a  uniform  biaxial  stress  system  N  and  T  at  infinity  (Fig.  8).  For  a  critical 
value  of  the  applied  loads  unstable  crack  extension  takes  place.  The  angle  of 
initial  crack  extension  ao  is  determined  by  assuming  that  the  crack  grows  in 
the  direction  of  the  maximum  circumferential  stress  cr#,  while  the  critical 


N  sT 


Fig.  8.  Geometry  of  a  rounded-off  angle  square  inclusion  partially  bonded  to  an  elastic 

matrix. 
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load  is  determined  by 

y^<Tjao)=K„  (52) 

where  K„  is  the  critical  stress  intensity  factor  which  is  a  material  constant. 

The  variation  of  the  crack  extension  angle  versus  the  crack  angle  /3  for 
s  =  -  I,  0  and  1  is  shown  in  Fig.  9,  while  Fig.  10  gives  the  dimensionless 
quantity  K„-KJ(TjJnR).  From  Fig.  10  the  critical  stress  T„  for  crack 
growth  is  determined. 


Fig.  9.  Variation  of  the  crack  extension  angle  a,  versus  half  crack  angle  p  for  a  square 
inclusion  for  s  =  - 1,  0  and  I. 


3S.  The  Triangular  Inclusion 

A  rigid  curvilinear  rounded-off  angle  triangular  inclusion  with  two 
locations  of  the  interfacial  crack  shown  in  Fig.  1 1  is  considered.  Using  the 
minimum  strain  energy  density  criterion  it  was  found  that  for  the  case  of 
Fig.  11(a)  the  crack  always  grows  from  its  tip  A.  The  variation  of  the 
dimensionless  quantity  =  A  \l(iSJ(Tl,R)  versus  half  crack  angle  w/2 
is  shown  in  Fig.  12.  Analogous  results  for  the  case  of  Fig.  1 1(b)  where  the 
crack  again  first  starts  from  its  tip  A  are  shown  in  Fig.  1 3.  From  these  figures 
the  critical  stress  for  unstable  crack  growth  is  determined. 
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Fio.  10.  Variation  of  the  dimensionless  stress  intensity  factor  versus  half  crack  angle  for 
a  squate  inclusion  for  s  =  —  1,  0  and  1. 


tN=sT 


tN=sT 


Fig.  II.  a  rounded-olT  angle  triangular  inclusion  partially  bonded  to  an  elastic  matrix. 
Geometrical  configuration  of  two  interfacial  crack  locations. 


4.  CONCLUDING  REMARKS 

The  failure  behavior  of  certain  particulate  composites  consisting  of  filler 
particles  with  elastic  moduli  much  higher  than  the  elastic  modulus  of  the 
matrix  was  studied.  Two  types  of  problems  modelling  the  composite  were 
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Fig.  12.  Variation  of  the  dimensionless  minimum  strain  energy  density  factor  versus  half 

crack  angle  (o/2  for  the  case  of  Fig.  11(a).  s=0,  0-25,  0-50,  0-75  and  10. 

considered:  (i)  rigid  inclusions  with  cuspidal  points  embedded  in  an  elastic 
matrix  and  (ii)  rigid  inclusions  partially  bonded  to  an  elastic  matrix.  In  the 
first  case  high  stress  concentrations  are  developed  in  the  vicinity  of  the 
cuspidal  points  which  constitute  nuclei  for  failure  initiation,  while  in  the 
second  case  failure  starts  from  the  tips  of  the  interfacial  crack  which 
coincides  with  the  unbonded  part  of  the  inclusion.  For  the  partially  bonded 
inclusion  a  general  solution  of  the  stress  and  displacement  fields  was 
obtained  for  any  inclusion  shape  using  the  method  of  complex  potentials. 
The  problem  was  reduced  to  a  Hilbert  problem  and  formulae  for 
determining  the  unknown  coefficients  of  the  solution  were  derived.  In  both 
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Fig.  13.  Variation  of  the  dimensionless  minimum  strain  energy  density  factor  versus  half 
crack  angle  (0/2  for  the  case  of  Fig.  11(b).  .v  =  0,  0  25.  0  50.  0-75  and  10. 


cases  results  were  obtained  for  special  inclusion  shapes  including  the  fiber, 
the  hypocycloidal,  the  astroidal,  the  square  and  the  triangular  inclusion. 
After  determining  the  stress  field  a  failure  analysis  of  the  composite  took 
place  using  the  maximum  circumferential  stress  criterion  and  the  strain 
energy  density  theory.  The  critical  load  for  fracture  initiation  from  the  more 
vulnerable  failure  sites  and  the  initial  fracture  angle  were  determined. 

The  results  of  this  work  shed  light  into  the  complicated  problem  of 
modelling  the  microstructure  of  particulate  composites  whose  reinforcing 
constituents  are  of  irregular  shape  like  the  various  inorganic  fillers,  the 
mefal  or  boron  filaments,  the  aggregate  or  sand  particles  in  concrete.  In 
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such  case  failure  of  the  composite  usually  initiates  from  the  sharp  angles  of 
the  inclusions  or  the  debonding  areas  of  the  different  phases.  The  analysis  of 
these  types  of  failure  mechanisms  is  of  major  importance  for  the  under¬ 
standing  of  the  failure  mode  of  the  composite. 
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ABSTRACT 

This  paper  derives  the  governing  equations  for  the  therw'-.mechanical 
behaviour  of  composites.  When  the  basic  equations  for  the  thermoelastic 
behaviour  of  solids  were  first  derived  in  the  nineteenth  century  several 
approximations  were  made.  The  effect  of  these  assumptions  are  discussed  and 
illustrated  by  the  results  of  a  simple  laboratory  test.  The  implications  of  this 
work  on  the  analysis  of  impact  damaged  laminates  are  then  discussed. 


I.  INTRODUCTION 

The  theory  describing  the  coupling  between  mechanical  deformation  and 
thermal  energy  of  an  elastic  body  was  first  published  in  1858  by  Lord 
Kelvin.  ‘  In  composite  materials  absorption  of  moisture  also  results  in 
internal  stresses  and/or  strains.  The  thermal  environment  may  also  interact 
with  moisture.  Indeed  it  is  generally  accepted  that  the  diffusion  of  moisture 
and  temperature  are  also  coupled.^'^  However,  tests  have  not  yet  been 
standardized  for  the  required  experimental  measurements.  The  current 
theories  for  composites  reduce  to  the  theory  of  classical  thermoelasticity 
when  moisture  effects  are  ignored  and  the  material  is  elastic.  However, 
when  deriving  the  expression  for  the  rate  of  change  of  entropy,  the  theory  of 
classical  thermoelasticity  assumes  that  the  stiffness  tensor  is  independent  of 
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temperature.  For  metals  this  assumption  has  been  shown  to  lead  to 
erroneous  results.* 

This  paper  presents  a  consistent  formulation  for  the  thermomechanical 
behaviour  of  composites  and  then  outlines  the  relevant  criteria  for  failure 
due  to  delamination  damage. 


2.  BASIC  EQUATIONS 

The  equations  governing  the  thermomechanical  behaviour  of  a  solid  body 
due  to  heating  and  external  forces  during  a  reversible  process  are  given 
below. 


2.1.  The  Constitutive  Equation 

Here  Cijki  is  the  stiffness  tensor,  jSy  and  (i>ij  are  coefficients  related  to  the 
thermal  and  moisture  expansion  coefficients  of  the  body  respectively  whilst 
7^  and  M,  are  the  reference  temperature  and  moisture  content  respectively. 


2.2.  Conservation  of  Mass  (continuity  equation) 

^  dx, 

"Df 


(2) 


Here  K  =  M  is  the  mass  flux  of  moisture  per  unit  volume  and  C  is  the 
‘instantaneous’  density. 


13.  Conservation  of  Momentum  (no  body  forces) 


(3) 


2.4.  Conservation  of  Energy 

E=a,jV,j-TS 

Here  S  is  the  entropy  and  £  is  the  internal  energy. 


(4) 
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2.5.  Rate  of  Change  of  Entropy 

TS=^-h,j  +  i:qM  (5) 

Here  q  is  the  heat  generated  when  1  g  of  moisture  is  absorbed  by  1  g  of 
the  material  and  h  is  the  heat  flux  tensor. 

With  this  formulation  the  specific  heat  C,  at  ifj  =  M  =  0,  is  defined  as 

~^u/cij.M  (6) 

The  remaining  equations  may  be  expressed  in  a  convenient  form  by 
introducing  the  free  energy,  F,  such  that 

F=E-TS  (7) 

The  quantity  F  can  be  expressed  in  terms  of  strain,  temperature  and 
moisture  in  a  manner  analogous  to  that  given  in  Ref.  5,  viz. 

F  =  l/2£y - p,je,j{T-  TJ -  (i>,^ £,/M -MJ  +  C^(T,M)  (8) 

Here  C,  is  a  function  of  temperature  and  moisture  and  as  shown  in  Ref.  5 


CC,= 


=  -r 


e^c^ 


(9) 


The  Duhamel-Neuman  law  tells  us  that 


whilst 


S  = 


dF 


(10) 


(11) 


With  this  notation  the  change  in  entropy  can  be  written  as 
.  dS.  dS  .  dS  . 


(12) 
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Now 


dT 


dT^ 


^CCJT 


and 


de,j 


whilst 


^  (Pij(T-  TJ)+  ^  {(t>,j{M-MJ) 

dS  ^  d^F 
dM  ~  dMdT  ~ 


(13) 

(14) 

(15) 


Substituting  eqns  (13),  (14)  and  (15)  into  (12)  we  finally  obtain  the 
‘Entropy  Equation’ 

rs=cc,t+4.r(^  rj)+ ^  (.^,,(m-mj) 

-  ^  (Cijk,)ek^  (16) 

If  we  now  substitute  eqn.  (5)  into  (16)  we  obtain 

-/i..(  +  CqM  =  CC,r 

(17) 

Fourier’s  law  for  heat  conduction  may  be  used  in  most  circumstances 
and  the  term  hj,j  occurring  in  eqn.  (17)  can  be  replaced  by 


where  k^J  is  the  thermal  conductivity  tensor. 

In  order  to  solve  the  above  set  of  equations  it  is  necessary  to  specify  an 
equation  describing  the  way  in  which  moisture  is  absorbed  into  a  com¬ 
posite.  Most  formulations  are  empirical.  The  equation  preferred  by  the 
authors  was  first  derived  in  Ref.  6  and  is  given  below,  viz. 

d  dM\  d 

where  Dy  and  X  are  material  constants.  An  extension  to  allow  for 
mechanical  coupling  in  the  matrix  material  replaced  the  term  Ton  the  right 
hand  side  of  eqn.  (18)  by  a  term  of  the  form  T+Ni]/  where  N  is  an 
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experimental  constant  and 

.A=ffu+AT-TJ  +  0(M-MJ  (19) 

Here  p  and  <f)  are  thermal  and  moisture  expansion  coefficients  for  the  matrix 
material. 


3.  IMPLICATIONS  FOR  CYCLIC  STRESSING 

If  an  elastic  body  is  subjected  to  cyclic  stresses  and  the  frequency  is  such  that 
the  process  is  adiabatic  then,  from  eqns  (5)  and  (16),  it  follows  that 

1=  -  ^Pij(T-TJ)--^{tf>,j{M-MJ)yj/i;C,  (20) 

Most  papers  neglect  the  term  dCuij/^T. 

If  this  is  done  is  said  to  be  related  to  T  by  a  material  constant. 
However,  it  is  well  documented’’®  that  even  for  metals  this  constant  is 
stress/strain  dependent. 

As  shown  in  Ref.  4  the  present  theory  is  able  to  accurately  predict  the 
stress  dependency  of  the  thermoelastic  constant  K  for  both  a  titanium  and 
an  aluminium  alloy.  The  thermoelastic  constant  and  the  Gruneisen 
parameter  for  a  metal  are  related  by  the  formulae 

/C=(l-2v)y/£  (21) 

In  order  to  illustrate  this  effect  let  us  consider  a  metal  bar  subjected  to 
a  uniaxial  stress 


(Ti,=S„  +  ASsincur,  (7(^  =  0  if  i,;Vl 

where  S„  is  the  mean  stress. 

Substituting  eqn.  (22)  into  (20)  gives 


(22) 


T  f  1  \  1  ^ 

CC,-  =  -  f  a-^  ^S,,ja>AScoswt+^  —  w(AS)^ sin 2(ut  (23) 

Integrating  with  respect  to  t  we  obtain 

AT  (  1  d£  „  \  A5  . 

=  —  I  a  —  sir 

T,  V  ^  /CCv 


sm  cut 


1  dE 


(AS)^  (1  —  cos  2cat) 


(24) 


where  7^  is  the  reference  temperature.  From  eqn.  (24)  we  see  that  the 
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thermoelastic  constant,  K,  is  given  by 

and  has  a  mean  stress  dependence.  Furthermore  if,  as  in  Ref.  4,  we  define 
K„=tx/CC,  then 


t  dK  _  \  dE 

K,W„~~Y^dT 


(26) 


The  quantity  (l/K<,)(5k/dS„)  was  measured  experimentally  by  Machin 
et  al.^  for  a  titanium  alloy  Ti-6A1-4V  and  an  aluminium  alloy  Al-2024  for 
which  the  values  dEldTwere  available.*'®  Table  1  lists  the  data  used  and  the 
comparison  between  the  theoretically  predicted  mean  stress  dependence 
and  the  experimental  results.  Good  agreement  between  theory  and 
experiment  is  clearly  evident. 

A  more  detailed  analysis  and  discussion  of  this  problem  is  contained  in 
Ref.  4. 


Table  1 

Comparison  of  theoretical  and  measured  mean  stress  dependence  of  K  (from  Ref  4 ) 


Material 

a 

rc-‘) 

£ 

(MPa) 

6E/dT 

(MPa/°C) 

(dK/dSJX;’ 

Theory 

(MPa-') 

Experiments® 

Ti-6A1-4V 

Al-2024 

90x  10'* 
2.3x10'’ 

o  o 

X  X 

-480 

-360 

4-33  X  10'* 
302x  10'* 

4-29  X  10* 
319x  10'* 

For  composite  materials  the  change  of  the  stiffness  tensor  with  temper¬ 
ature  ‘  “  is  far  greater  than  for  metals,  see  Table  2,  and  should  not  be  ignored. 


Table  2 

Graphite-epoxy  ( ASI350I-6)  orthotropic  elastic  properties  (from  Ref.  10) 

Temperature 

(T) 

Parameter  (lb/in‘ 

xlO'®) 

Vu 

£T 

Ej  and  Ej 

G.2 

—  65  (dry) 

0-3 

18-70 

17-60 

2-30 

1-00 

RT  (dry) 

0-3 

18-70 

17-60 

1-90 

0-85 

220  (dry) 

0-3 

18-40 

17-30 

1-65 

0-65 

250  (dry) 

0-3 

18-40 

17-20 

1-60 

0-63 

220  (wet) 

0-3 

17-30 

1-32 

0-52 
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The  temperature  rise  due  to  mechanical  loading  is  usually  small  (i.e.  of 
the  order  1/K)  unless  there  is  extensive  plasticity.  However  in  the  vicinity  of 
a  crack  or  a  delamination  c^j  is  infinite  and  eqn  (16)  predicts  that  the  glass 
transition  temperature  will  be  exceeded  in  a  small  localized  region. 
However  the  validity  of  this  approach  which  is  based  on  reversible 
thermodynamics  is  questionable  once  crack  extension  or  delamination 
growth  has  occurred. 

The  residual  thermal  strains  ej-  due  to  the  curing  process  may  also  play 
a  role  in  the  variability  of  fatigue  Iife.The  residual  strains  may  be  interpreted 
as  the  mean  strain  level  about  which  the  mechanical  strains  oscillate.  As 
shown  above  the  temperature  depends  both  on  the  mean  stress/strain  and 
the  cyclic  stress/strain  fields. 


4.  FRACTURE  MECHANICS 

There  are  two  fracture  parameters  which  are  widely  used  to  predict  the 
residual  strength  of  delaminated  composite  laminates.  These  are; 

(i)  the  energy  release  rate  approach; 

(ii)  the  strain  energy  density  approach. 

4.1.  The  Energy  Release  Rate  Approach 

For  mode  1  self  similar  crack  growth  of  a  through  crack  in  the  absence  of 
body  forces  the  energy  release  rate  G  can  be  written  as 

C  =  £(»'•», (27) 

where  W  is  the  energy  density,  F,  is  a  vanishing  small  closed  path  around 
the  tip  with  normal  n  and  are  the  components  of  the  traction  vector 
on  the  path. 

Here  W  is  defined  as 

W=  l/2e„Qj„e,j-p,je,j{T-  +  M)  (28) 

Let  us  now  consider  the  integral  J,  which  we  will  define  as 


where  F,  is  the  external  boundary  of  the  body.  As  mentioned  in  Ref.  11 
there  is  a  tendency  to  drop  the  subscript  s  and  refer  to  J,  as  J.  Using 


56 


R.  Jones.  T.  E.  Toy  and  J.  F.  Williams 


Green’s  theorem  it  follows  that 


(30) 


where  —  V\  is  the  area  between  the  curves  F,  and  F,.  If  we  assume  that  Pij, 
(pij  and  Cij^,  are  only  functions  of  T  and  M  then  the  term  dvi/jdx^  can  be 
written  as 


and 


dW 

dx 


,  ^x,V  + 


+  <^y(M-MJ)+|^jdK 


(31) 


(32) 


Thus  y,  will  not  equal  the  energy  release  rate  G  unless  the  area  integral 
vanishes.  At  constant  moisture  and  temperature  J,  may  be  equal  to  G. 
HoweVer  in  general  the  area  integral  will  be  non-zero  and  J^,  which  is 
measured  experimentally  from  the  movement  of  the  load  points, will  not 
equal  G. 

In  service  aircraft  heating  is  often  localized  and  the  moisture  content 
varies.  This  will  produce  a  spatial  variation  in  the  tensor  Cy*/  with  the  result 
that  the  area  integral  will  in  general  be  non-zero.  Consequently  when 
designing  laboratory  tests  care  should  be  taken  to  reproduce  the  near  tip 
stress,  strain,  temperature  and  moisture  fields  rather  than  reproducing  the 
‘global  behaviour’.  This  is  particularly  true  if  the  aim  of  the  test  is  to 
establish  such  quantities  as  the  critical  damage  size  or  the  maximum 
permissible  load. 

For  a  three-dimensional  fracture  problem  the  integral  on  the  right  hand 
side  of  eqn.  (30)  is  no  longer  equal  to  the  energy  release  rate  and  is  referred  to 
as  T*. 
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5.  STRAIN  ENERGY  DENSITY 

In  the  strain  energy  density  approach  failure  is  assumed  to  occur  when  the 
available  energy  density  W' ,  at  a  distance  in  front  of  the  delamination  in 
the  direction  of  growth  reaches  a  critical  value  W^.  The  value  is 
dependent  both  on  the  values  of  dK(=£j,  +  633)  and  d/4,  the  area  of 

crack  growth. 

For  thermomechanical  problems 

=  l/2a,,e,^  -  TJ  -  -  M„)  -  W,  (33) 

where  =  <l>ij  =  ilfijCijia  and  Wf  is  the  energy  density  in  the  fibre. 

As  an  illustration  of  this  approach  consider  an  impact  damaged  laminate 
with  a  fastener  hole  under  compression.  The  dimensions  of  the  model  are 
the  same  as  those  used  in  the  experimental  work  of  Ref  1 3,  see  Fig.  1  in  this 
reference.  The  specimen  tested  was  a  [O/45/O2/-45/O2/— 45/0],  T300/5208 
graphite-epoxy  laminate  and  contained  a  centrally  located  hole  9-5  mm  in 
diameter,  surrounded  by  delamination  damage  due  to  impact  and  poor 
drilling.  The  elements  used  are  mostly  twenty-noded  isoparametric  ele¬ 
ments  with  directionally  reduced  integration  and  2x2x3  Gaussian 
quadrature  points,  with  the  3  points  being  taken  through  the  ply  thickness. 
The  crack  tip  elements  along  the  circular  delamination  are  fifteen-noded 
isoparametric  wedge  elements. 

The  initial  damage  around  the  fastener  hole  (from  Ref  13)  is  modelled  as 
a  circular  delamination  of  radius  1 3-75  mm  between  the  second  and  third 
plies  (i.e.  between  45°  and  0°  plies).  It  can  be  seen  from  the  ultrasonic  C-scan 
that  the  initial  delamination  is  nearly  circular.'^ 

The  two  plies  above  and  below  the  delamination  and  the  matrix  region 
around  the  delamination  are  modelled  separately  with  ordinary  three- 
dimensional  elements  while  the  remaining  20  plies  are  modelled  with 
super-elements  with  displacements  varying  quadratically  in  the  local 
isoparametric  co-ordinate  system.  The  material  properties  used  are  those  of 
ASl/3501-6. 

It  is  important  that  in  the  FE  model,  the  faces  of  the  delamination  are 
prevented  from  overlapping.  Otherwise,  non-physical  solutions  may  be 
obtained.  By  examining  the  solutions  of  the  displacements,  it  is  found  that 
some  parts  of  the  delaminated  faces  have  overlapped.  Thus  a  series  of 
constraint  equations  are  applied  to  appropriate  nodes  to  simulate  local 
closure. 

Two  load  cases  are  considered,  viz. 
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1.  A  compressive  load  of  150  kN  applied  at  the  ends  of  the  model. 

2.  The  load  case  mentioned  above  together  with  a  10%  moisture 
content.  To  induce  stresses  due  to  moisture  the  external  boundary  of 
the  specimen  is  prevented  from  moving  in  the  direction  perpendicular 
to  the  primary  load  (i.e.  load  case  1).  Table  3  shows  the  maximum 
value  of  T*  for  the  two  load  cases  as  well  as  the  maximum  values  of 
the  parameter  W^  where 

IT,  =  limit  ITn.ds  (34) 

Jr, 


Table  3 

Maximum  values  of  T*  and  W,(  J/M) 
around  the  delamination 


Load  case 

T* 

IT, 

1 

71-8 

541 

2 

1680 

51-3 

Here  we  see  that  the  value  of  T*  is  increased  dramatically  by  the  presence 
of  moisture.  However,  for  the  present  problem  the  interlaminar  stresses  or., 
T„  and  Tyj  are  relatively  unaffected  by  the  moisture  content.  This 
phenomenon  is  also  seen  in  the  ratio  of  the  maximum  value  of  the  strain 
energy  density  IT  in  the  matrix  material  directly  in  front  of  the  delamina¬ 
tion  for  the  two  load  cases: 

IT(load  case  l)/IT(load  case  2)=  102 

Indeed  from  Table  3  we  see  that  W,,  which  involves  the  integral  of  the 
energy  density  around  the  delamination,  is  also  relatively  unaffected  by  the 
presence  of  moisture.  This  infers  that,  in  the  present  problem,  the  presence 
of  moisture  does  not  significantly  increase  the  likelihood  of  static  failure  by 
means  of  delamination  growth  along  the  interface.  This  observation  is 
consistent  with  the  experimental  results  given  in  Refs.  14  and  1 5.  We  also  see 
that  if  the  growth  of  the  damage  is  likely  to  be  non-self-similar  then  energy 
release  rate  methods  must  be  used  with  extreme  caution. 
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ABSTRACT 

A  method  of  calculating  the  energy  release  rate,  G,for  plane  delaminations  in 
composites  is  given.  This  is  couched  in  terms  of  local  moments  and forces  and  is 
a  useful  general  method  of  obtaining  such  solutions.  A  scheme  for  the  exact 
partitioning  of  the  G  values  into  mode  I  and  II  components  is  also  given.  The 
method  is  then  used  to  give  results  for  both  a  constant  ratio  mixed  mode  test 
and  one  in  which  the  ratio  continuously  varies.  Data  are  presented  for  three 
matrix  resins  used  in  carbon  fibre  laminates  and  for  two  of  them,  epoxy  and 
PEEK,  a  unique  locus  of  G,  and  G„  is  defined.  For  bismaleimide  fibre 
bridging  results  in  G^  increasing  with  crack  growth  and  consequent  differences 
between  the  tests.  The  test  methods  are  generally  confirmed  as  satisfactory. 


1.  INTRODUCTION 

Laminates  of  high  strength  and  stiffness  fibres  with  matrices  of  tough 
polymers  are  of  considerable  commercial  interest  since  they  offer  a  sig¬ 
nificant  potential  weight  advantage  over  conventional  materials.  Their  very 
nature  results  in  a  structure  in  which  these  properties  are  realised  in  the 
plane  of  the  laminate  but  their  weakness  lies  in  the  through  thickness  or 
translaminar  direction.  Here  strength  and  stiffness  are  approximately  those 
of  the  matrix.  In  designing  structures  care  is  taken  to  exploit  the  in-plane 
properties  but  inevitably  some  forms  of  loading  can  induce  failure  between 
laminate  layers,  i.e.  delamination,  and  this  is  often  the  limiting  property  of 
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the  structure.  It  is  therefore  important  to  have  a  good  understanding  of 
delamination  and  to  be  able  to  characterise  it  under  any  loading  system  so 
that  design  predictions  can  be  made.  This  paper  describes  a  scheme  by 
which  such  a  characterisation  may  be  carried  out. 

The  methodology  is  essentially  that  of  linear  elastic  fracture  mechanics 
(LEFM)  in  that  the  structure  behaves  in  a  linear  elastic  fashion  and  failure  is 
assumed  to  occur  in  a  plane  parallel  to  the  surface.  This  is  a  very  important 
simplification  in  the  analysis  since  it  avoids  the  necessity  of  predicting  the 
direction  of  crack  growth  under  complicated  loading  systems.  It  is  assumed 
that  the  laminated  structure  is  such  that  the  crack  propagation  is  forced  to 
occur  in  the  plane  of  the  sheet.  The  loadings  which  effect  this  growth  can  be 
bending  moments,  shear  forces  or  in-plane  loads  and  may  involve  buckling 
of  the  sheet.  In  such  cases  the  crack  can  be  forced  to  grow  under  a 
combination  of  an  opening  mode  (mode  1),  and  a  sliding  mode  (mode  II).  In 
isotropic  homogeneous  systems  such  loading  can  lead  to  non-colinear 
growth  but  here  it  is  forced  to  be  colinear.  Under  these  circumstances  we 
must  seek  a  fracture  criterion  for  mixed-mode  crack  propagation  and  to  do 
this  we  must  have  tests  which  enable  us  to  propagate  cracks  under 
mixed-mode  conditions.  In  addition  we  must  be  able  to  analyse  the  loading 
condition  so  that  we  may  determine  the  energy  release  rates  for  each  mode. 
Such  an  analysis  coupled  with  mixed-mode  fracture  criteria  will  provide 
a  basis  for  a  design  method. 


2.  METHOD  OF  ANALYSIS 

In  this  analysis  we  shall  derive  expressions  for  the  energy  release  rate  for  the 
situation  shown  in  Fig.  1  in  which  there  is  a  single,  through  thickness, 
delamination  in  a  laminate  of  thickness  2h  which  is  located  a  distance  h^ 
from  one  surface  and  in  which  hi  +  h2-2h,  h^<h2-  This  is  essentially 
a  one-dimensional  model  and  enables  the  analysis  to  be  conducted  in  terms 


Fig.  I.  Delamination. 
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of  simple  beam  theory.  Practical  problems  of  other  shapes  could  be  tackled 
in  a  similar  way  but  would  lead  to  more  complicated  solutions.  The  general 
method  is  dealt  with  in  detail  elsewhere and  will  be  considered  only  in 
outline  here.  Let  us  consider  one  end  of  the  delamination  of  length  a  as 
shown  in  Fig.  2.  The  general  loadings  at  this  point  are  a  moment  of 


Fig.  2.  Loading. 


+  M2,  a  shear  force  of  Q1  +  Q2  and  an  axial  load  of  Pi  H-P^.  In  the 
cracked  section  these  loads  are  A/,,  6,  and  P,  on  the  upper  section  of 
thickness  h^  and  similarly  Mj,  ^2  on  /i2-  Let  us  consider  first  only  the 
moments  and  note  that  for  a  section  of  second  moment  of  area  /  and  axial 
modulus  £  the  strain  energy  per  unit  length  of  the  beam  is  given  by 
In  Fig.  3  we  have  the  moment  versus  angle  of  rotation,  4>,  relationship  for 


Fig.  3.  Bending  moment  loading  lines. 


the  section  with  a  crack  of  length  a,  i.e.  the  line  OA,  which  is  linear  in  this 
case.  If  the  crack  grows  to  (o  +  da)  then  the  loading  line  becomes  OA'  and 
the  change  in  energy  of  the  system  for  this  change  in  a  is  given  by  the  shaded 
triangular  area  OAA'. 

Now  the  strain  energy  l/,=  \/2{M(j))  so  G  may  be  written  as: 

1  (  dM\  1 

dU  =  -M(p  +  {M  +  —m--{M  +  dM)((t)  +  d<t>) 
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i.e. 


or 


BG 


_dU_l/^d<f> 
da  2  \  da 


^  da) 


1  .dM 

BG= 

2  da 


$  const 


=  + 


A#  const 


da 

dU, 

da 


^  const 


M  const 


(1) 


The  former  is  the  more  usual  route  but  the  latter  is  more  convenient  here.  It 
should  be  noted  that  the  use  of  constant  or  M  in  the  analysis  to  find 
G  does  not  imply  that  this  is  the  loading  in  practice  but  is  simply  a  device  for 
calculation.  Returning  now  to  Fig.  2  we  can  compute  dUJda  at  constant 
M  here  by  calculating  the  energy  change  within  the  contour  ABB'A  when 
a  increases  by  da  (O  to  O'); 

dU,  \  M\  1  Mi  1(Mi  +  M2)^ 
da~2  Ell'll  EI2  2  £/o 

dl3 

where  /,  =  8(1 /o  =  8/ 

,  Bh^ 

and  ^  =  ‘^12 


The  expression  for  the  total  G  now  becomes, 


16B£/ 


Mi  Ml 


-(M,+M/ 


(2) 


This  is  an  important  result  since  it  enables  G  to  be  calculated  from  the  local 
values  of  M,  B,  E  and  /  of  the  crack  tip  without  recourse  to  the  values  in 
remote  regions.  Similar  analyses  lead  to  expressions  for  G  for  P  and  Q; 


(3) 

and 

(4) 

when  A  =  Bh  and  n  is  the  shear  modulus. 

These  expressions  are  for  the  total  energy  release  rate  and  it  is  necessary 
to  partition  them  into  modes  I  and  II.  In  bending,  mode  II  occurs  when 
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there  is  crack  growth  with  the  curvature  in  both  sections  equal.  If  we  have 
a  moment  Mn  on  /i,  and  ij/Mn  on  hj  then  we  have  pure  mode  II  when: 

d(f>  Mit  ij/Mii 

da”E/7^  EI2 
i.e. 


Mode  I  requires  equal  moments  in  opposite  senses  so  we  may  write; 

Afi  =  Mii  — iW,  and  M2  =  i/'Mii  +  M| 

Substituting  these  expressions  in  eqn.  (2)  we  have  a  term  in  Mu,  one  in 
Ml  but  that  in  A/|Mn  is  zero  thus  giving  exact  partitioning.  We  may  thus 
write; 

Af,\  (1+iA)  1 

'  BEI  ’  BEI  ' 

r  +  3  (l-j) 

"  BEI  16  ^  BEI  + 

Axial  loads  give  only  mode  11  so  that;  G,  =  0  and  G„  is  given  by  eqn.  (3)  and 
for  shear  forces  G„  =  0  and  G,  is  given  by  eqn.  (4). 


3.  MIXED  MODE  TESTS 

The  analysis  of  mixed  mode  crack  propagation  requires  some  form  of 
criterion  if  it  is  to  be  employed  in  design.  The  most  obvious  would  be  that  the 
total  energy  release  rate  remained  the  same,  i.e.  G  constant,  but  this  appears 
not  to  be  so  and  there  is  usually  a  distinct  difference  between  G,c  and  G,ic. 
The  criterion  can  be  represented  as  a  locus  of  G,  versus  G„  at  fracture  and 
some  form  of  relationship, 

F{Gi,  G„)  =  0 

is  required.  In  order  to  explore  this  it  is  necessary  to  perform  tests  in  various 
combinations  of  modes  I  and  II  and  the  following  have  been  employed: 

(1)  pure  mode  I; 

(2)  pure  mode  II; 
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(3)  mixed  mode  I/II  in  which  the  ratio  Gi/G,,  remained  constant; 

(4)  mixed  mode  I/II  in  which  G^/G„  varied  continuously. 

These  latter  two  enable  the  important  characteristic  of  history  dependence 
to  be  studied. 

Figure  4  shows  the  test  configurations  used  in  this  test  series  and  the 


PI/L 

3)  Fixed  Ratio  Mixed  Mode  4)  Variable  Ratio  Mixed  Mode 

Fig.  4.  Test  configurations. 


analysis  employed  is  as  follows: 


(1) 

Af2  =  M,= 

Pa,  = 

^2  '  BEl 

G„  =  0 

(6) 

(2) 

Pa  ,  1 

r  _  3 

16  BE/ 

(7) 
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(3)  M,=:0,  M2  =  Pa 

_  1 
G,= 


BEI  16(l-^)^(l  +  iA)’ 


G,  3(1-^)'^ 
G„ 


,  constant 


G„  = 


PV  3(1 -^) 
BEI  I6^\l+il/) 


(8) 


4.  EXPERIMENTAL  RESULTS 

The  experiments  will  be  described  in  detail  elsewhere  ^  *  and  only  an  outline 
of  the  results  will  be  given  here.  Ail  four  test  configurations  were  used  on 
unidirectional  laminates  of  carbon  fibres  using  three  different  matrices,  i.e. 
poly(ether-ether  ketone)  (PEEK),  an  epoxy  and  bismaleimide.  In  each  test 
the  cracks  were  grown  by  loading  Che  specimens  in  an  Instron  testing 
machine  and  recording  the  load,  load  point  displacement,  and  the  crack 
length  by  observing  marks  on  the  specimen  edge.  G  values  can,  of  course,  be 
found  via  eqns  (6-9)  from  P  and  a  only  but  the  deflection,  S,  is  useful  since 
for  the  mode  I  test  with  a  uniform  section  El  can  be  found  from  the 
relationship; 

.  2Pa^ 

and  this  value  was  used  in  all  the  subsequent  tests. 

During  each  test  it  was  possible  to  evaluate  G  as  the  crack  grew  and  Fig. 
5  shows  G  as  a  function  of  crack  growth  for  different  modes  of  testing  for 
PEEK.  In  mode  I  there  is  clearly  a  constant  G  value  but  in  mode  II  there  is 
an  increase  in  G  as  u  increases.  There  is  an  extreme  example  of  this  effect 
shown  in  Fig.  6  for  bismaleimide  and  this  gives  some  hint  as  to  the  source 
since  adhesion  of  the  fibres  to  the  matrix  here  was  not  good  (the  fibres  were 
not  correctly  primed).  This  results  in  the  fracture  occurring  in  several  planes 
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Fig.  5.  Variation  of  G  with  crack  growth  for  PEEK 


and  in  fibres  ‘bridging’  the  two  arms  giving  higher  values  which  increase  as 
the  crack  grows.  In  the  PEEK  and  epoxy  specimens  the  effect  is  much  less 
marked  because  the  adhesion  is  better. 

Results  for  the  three  materials  as  C,  versus  G„  are  shown  in  Figs  7, 8  and 
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G„  (kj/m2) 


Fig.  7.  Failure  locus  under  mixed  mode  loading  for  PEEK. 


Gii  (kj/m^) 

Fig.  8.  Failure  locus  under  mixed  mode  loading  for  epoxy. 
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FiC.  9.  Failure  locus  under  mixed  mode  loading  for  bismaleimide. 

9.  For  PEEK  and  epoxy  there  is  a  quite  well  defined  locus  for  both  forms  of 
mixed  mode  test  and  there  is  little  evidence  of  fibre  bridging.  For  the 
varying  ratio  tests  l/L  was  varied  to  check  independence  of  test  method  and 
the  correctness  of  the  analysis.  The  agreement  between  the  two  forms  of  test 
is  encouraging.  For  bismaleimide  the  fibre  bridging  results  in  large 
discrepancies  and  the  notion  of  a  unique  locus  cannot  be  used.  Since  one 
would  hope  for  good  matrix-fibre  adhesion  in  practice  it  is  felt  that  the 
single  line  is  likely  to  prove  useful.  The  form  of  the  curve  is  somewhat  below 
a  linear  relationship  between  G,c  and  G„c- 
The  quality  of  the  data  confirms  that  the  test  configurations  and  the 
analysis  employed  are  satisfactory  and  represent  a  useful  approach  for 
analysing  this  form  of  fracture. 
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ABSTRACT 

A  procedure  has  been  developed  which  enables  one  to  employ  a  computer 
program  to  carry  out  the  essential  computations  required  for  the  generation  of 
non-linear  constitutive  expressions  for  the  cases  where  the  material  con¬ 
sidered  belongs  to  one  of  the  32  crystal  classes.  This  procedure  is  outlined  here 
and  extended  so  as  to  also  cover  the  cases  where  the  group  defining  the 
material  symmetry  is  one  of  the  five  transverse  isotropy  groups. 


1.  INTRODUCTION 

Composite  materials  are  replacing  more  traditional  materials  in  a  wide 
variety  of  applications.  The  objective  is  to  design  a  material  which  has  the 
properties  required  for  the  given  application.  We  may  take  a  first  step  in  this 
direction  if  we  can  specify  the  form  of  the  constitutive  equation  which 
would  best  suit  out  purposes.  The  form  of  a  constitutive  equation 
describing  the  response  of  a  material  is  essentially  dictated  by  the  symmetry 
properties  of  the  material.  We  may  search  through  the  various  categories  of 
constitutive  equations  associated  with  materials  which  possess  symmetry 
properties  in  order  to  find  which  most  closely  matches  the  desired  form.  We 
may  then  hope  to  fabricate  a  material  with  the  appropriate  symmetry  by 
judicious  placement  of  fibers  in  a  matrix  comprised  of  an  isotropic  material. 
In  order  to  aid  in  the  recognition  process,  we  write  the  constitutive 
expression  in  matrix  form  where  the  matrices  describing  the  material 
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properties  are  diagonalized  to  the  extent  possible.  We  employ  methods 
from  group  representation  theory  to  attain  this  end.  We  plan  to  consider 
non-linear  constitutive  expressions  which  generally  require  a  substantial 
amount  of  computation.  To  overcome  any  difficulties  in  this  direction,  we 
have  developed  a  computer  program  which  will  generate  the  diagonalized 
matrix  form  of  a  wide  class  of  constitutive  expressions  for  materials 
possessing  the  symmetry  associated  with  any  of  the  crystallographic 
groups.*'^  In  this  study,  we  give  the  basic  information  required  to  extend 
these  results  to  the  cases  where  the  material  possesses  symmetry  properties 
characterized  by  one  of  the  live  transverse  isotropy  groups  which  we  denote 
by  r„  Tj, . . . ,  Tj.  The  group  T,  defines  the  symmetry  of  a  material  which 
possesses  rotational  symmetry  about  the  Xj  axis.  The  group  T^  defines  this 
symmetry  of  a  material  which  possesses  rotational  symmetry  about  the  Xj 
axis  and  for  which  the  plane  containing  the  Xj  and  Xj  axes  is  a  plane  of 
symmetry.  The  group  Tj  defines  the  symmetry  of  a  material  possessing 
rotational  symmetry  about  the  Xj  axis  and  for  which  the  plane  containing 
the  Xj  and  Xj  axes  is  a  plane  of  symmetry.  The  group  T*  is  associated  with 
a  material  possessing  rotational  symmetry  about  the  Xj  axis  and  for  which 
the  planes  containing  the  Xj  and  Xj  axes  and  the  x,  and  Xj  axes  are  both 
planes  of  symmetry.  The  group  Tj  is  associated  with  a  material  possessing 
rotational  symmetry  about  the  Xj  axis  and  for  which  the  Xj  axis  is 
a  two-fold  axis  of  symmetry. 

We  observe  that  the  consideration  of  transversely  isotropic  materials 
may  arise  in  connection  with  composites  in  the  following  ways.  Consider 
a  circular  cylinder  in  which  is  embedded  a  number  of  fibers  which  are 
parallel  to  the  axis  of  the  cylinder.  If  the  fibers  are  uniformly  distributed 
over  the  cross-section  of  the  cylinder,  we  may  assume  that  the  material 
possesses  rotational  symmetry  about  the  axis  of  the  cylinder.  In  addition, 
the  plane  perpendicular  to  the  axis  of  the  cylinder  and  the  planes  containing 
the  axis  of  the  cylinder  are  planes  of  symmetry.  The  group  defining  the 
material  symmetry  would  be  the  transverse  isotropy  group  denoted  by  T*. 
One  may  fabricate  a  material  consisting  of  layers  of  thin  sheets  in  which  an 
array  of  parallel  fibers  is  present.  If  the  fibers  in  any  given  layer  make  an 
angle  of  2n/n  with  the  fibers  in  the  adjacent  layers,  we  may  consider  the 
composite  to  have  an  n-fold  axis  of  symmetry  perpendicular  to  the  layers. 
Further  we  may  assume  that  the  material  possesses  a  plane  of  symmetry 
which  contains  the  n-fold  axis  of  symmetry.  If  n  becomes  large,  this  would 
render  the  material  approximately  transversely  isotropic.  The  group 
defining  the  material  symmetry  would  be  the  transverse  isotropy  group  Tj. 

We  discuss  below  the  procedure  employed  to  generate  the  block  diagonal 
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form  of  the  matrix  constitutive  equations.  We  then  give  the  information 
required  to  employ  this  procedure  for  the  transversely  isotropic  groups  Tj, 
Tj, . . .  .Tj.  We  finally  give  examples  of  the  application  of  these  results  to  the 
generation  of  non-linear  constitutive  expressions.  The  work  discussed 
below  constitutes  part  of  the  doctoral  dissertation^  of  one  of  the  authors 
(G.  Bao). 


2.  TRANSVERSELY  ISOTROPIC  MATERIALS 

The  constitutive  expressions  which  we  consider  are  tensor-valued  functions 
of  one  or  more  tensors  Sj,  Sj  , . . .  of  degrees  n,,  , . . .  in  these  tensors 

which  are  invariant  under  a  group  F  defining  the  material  symmetry.  We 
are  primarily  interested  in  the  cases  where  the  material  symmetry  is  defined 
by  one  of  the  five  groups  . . . ,  T,  associated  with  the  various  types  of 
transversely  isotropic  materials.  These  groups  are  defined  by  specifying  the 
groups  of  3  X  3  matrices  which  define  the  set  of  equivalent  reference  frames 
associated  with  the  material.  Thus  the  group  T^  is  comprised  of  the  matrices 


cos  9 

sin0. 

0 

—  sin  9, 

COS0, 

0 

,  0^0 <271 

(1) 

0, 

0, 

1 

.Let  e,  denote  the  base  vectors  associated  with  rectangular  Cartesian 
coordinate  system  x.  Let  Cj  denote  the  base  vectors  associated  with  the 
reference  frame  x  which  is  obtained  by  rotating  the  reference  frame 
X  through  9  radians  counter-clockwise.  We  have 

=  (2) 

where  the  matrix  Q(0)=  Il6j/0)ll  is  defined  by  (1).  If  x  and  x  are  equivalent 
reference  frames,  i.e.  if  x  arises  from  x  by  applying  a  symmetry  operation  to 
X,  the  constitutive  equation  is  required  to  have  the  same  form  when  referred 
to  either  reference  frame.  Let  and  Ejj  denote  (absolute)  second-order 
tensors.  Then,  if 

~  ■  ■  ■  ^mn  (3) 

is  the  constitutive  expression  when  referred  to  the  x  frame,  the  constitutive 
expression  where  referred  to  the  x  frame  is  given  by 

T  =C  F  F 


(4) 
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where 


T^iJ  QipQjqTpq> 
^kt  —  QksQlt^si 


^ij..,,=  QipQjq  •  •  Qnr^pq.r 


(5) 


If  the  reference  frame  x  and  x  are  equivalent,  we  require  that 


=  C 

'ijkl...mn  ^  ijkl...mn 


(6) 


If  the  symmetry  operations  consist  of  all  rotations  about  the  Xj  axis,  the 
property  tensor  „„  must  satisfy  the  equations 


Qipmjqm-Qnuie)Cpq....=c,j.,,  (7) 

for  0  <  0  <  2n.  We  say  that  the  tensor  C,j  „  which  satisfies  eqn  (7)  for  0  ^  0  < 
2n  is  invariant  under  the  group  Tp  We  may  express  a  tensor  which  is 
invariant  under  the  group  T,  as  a  linear  combination  of  the  outer  products 
of  the  fundamental  tensors 


^3i'  +  — “iV'  ^  U^2j~  Pij  (8) 

Thus,  one  may  express  the  general  fourth-order  tensor  which  is  invariant 
under  the  group  T,  as  a  linear  combination  of  the 

^3i^37*lkl*^'  (^) 

PijPkl’ 

where  the  numbers  following  the  tensors  denote  the  number  of  distinct 
isomers  of  the  tensor.  An  isomer  of  a  tensor  (7,^^,  is  obtained  by  permuting 
the  subscripts  i,  j,  k,  I  of  the  tensor.  We  have  noted  that  a.j  =  oij^  and 
=  —  Pjf.  We  further  note  that 

^iA/  =  “i/a;*  (10) 

and  that  only  three  of  the  six  isomers  of  are  linearly  independent. 
Thus,  we  have 

“.A/  +  “i*4  +  “A  =  0 

«iA<  +  ®*A  +  “i/i*=0  (11) 

ai»)5;/  +  “;*^/i  +  “/»^O  =  0 

The  existence  of  relations  such  as  (1 1)  renders  the  generation  of  the  general 
tensors  of  orders  5, 6, . . .  which  are  invariant  under  T,  a  non-trivial  matter. 
This  may  be  accomplished  using  a  method**  which  employs  Young 
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tableaux.  However  the  procedure  for  generating  the  form  of  a  constitutive 
equation  based  on  listing  the  general  form  of  the  property  tensor  Cy*,  „„ 
and  then  substituting  into  (3)  would  generally  prove  to  be  cumbersome. 

It  is  preferable  to  employ  a  procedure  based  on  group  representation 
theory.  A  set  of  matrices  P(0)  which  is  in  one  to  one  correspondence  with  the 
matrices  Q(fl)  comprising  the  group  Tj  and  such  that  P(0i)P(02)  cor¬ 
responds  to  Q(0i)Q(02)  is  said  to  form  a  matrix  representation  of  the  group 
Tp  The  set  of  matrices  KP(0)K"‘  where  det  K#0  also  forms  a  matrix 
representation  of  Tj  which  is  said  to  be  equivalent  to  the  representation 
P(0).  An  appropriate  choice  of  the  matrix  K  enables  us  to  write 

KP(0)K-*  =  a,P,(0)-ha2P2(0)-h  •  •  (12) 

in  block  diagonal  form  where  aj,  a2, . . .  are  positive  integers  and  where 

P,(0),  0,  0 

2P,(e)  +  P^i0)=  0,  P,(0),  0  (13) 

0,  0,  P2(0) 

We  say  that  the  representation  P(0)  may  be  decomposed  into  the  direct  sum 
of  the  representations  P,  (0),  P2(0), ....  If  a  representation  P(0)  cannot  be 
decomposed,  it  is  referred  to  as  an  irreducible  representation.  The 
irreducible  representations  associated  with  T,  are  all  one-dimensional  and 
are  defined  by  listing  the  1  x  1  matrix  corresponding  to  the  matrix  Q(0).  We 
define  these  representations  below. 

Vo:  » 

Vp:  (p=l,2,...)  (14) 

Tp;  e‘P®  (p=l,2,...) 

Vq  denotes  the  identity  representation  where  the  same  number  1  cor¬ 
responds  to  each  Q(0);  denote  the  representation  where  e"'*^  corresponds 
toQ(0),...  . 

We  consider  the  manner  in  which  a  vector  transforms  under  the  group 
Ti .  The  components  X;  of  a  vector  x  when  referred  to  the  reference  frame 
X  with  base  vectors  Cj  =  Qij(d)ej  are  related  to  the  components  x,  of  x  when 
referred  to  the  x  frame  with  base  vectors  e,  by  the  equations 

X,  cos  0,  sin  0,  0 

Xj  =  Qij{0)Xj  or  X2  =  —  sin  0,  cos  0,  0 
Xj  0,  0,  I 


(15) 
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With  (15)  we  readily  see  that 

Xj+iXj  e“‘*,  0,  0  Xi+iXj 

Xj— ixj  =  0,  e'*,  0  X,  — ixj  (16) 

X3  0,  0,  1  Xj 

This  tells  us  that  the  transformation  properties  of  Xj  +ix2,x,  —  ixj,  Xj  under 
the  group  Tj  are  defined  respectively  by  the  irreducible  representations  y,, 
Tj  and  respectively.  We  immediately  see  that  the  transformation 
properties  of 


(Xi  +iX2)^  (X,  +iX2K.Xi  -iXj),  (x,  +iX2)X3, 

(Xi-iXj)^  (X,-iX2)X3,  xl 


(17) 


are  defined  by  the  irreducible  representations  yj*  Vo*  Vi-  Tj,  Fi,  '/„ 
respectively.  Since  the  components  x,Xj  transform  in  the  same  manner  as  do 
the  components  of  a  symmetric  second-order  tensor,  we  see  that  the 
transformation  properties  of 


^11  ^11  +^22*  ^13  +  >^23' 

■^11  ~^22“'2iSi2,  5,3  — iS23,  ^33 


(18) 


transform  according  to  y2,  y’o,  >’i,  r2,  F,,  y^  respectively.  Similarly  we  see 
that 


(s,  1  —  522+215,2)^,  (s,,  —S22 +215,2X5, 1  +522), 

(5, ,-522  +  215,2X5, 3  +  1523),  (19) 

(5, ,  —  522  +  215,2X5, ,  —  522  ~  215, 2),  •  •  • 

transform  according  to  y^,  y2,  >’3,  ’/q,  . . .  respectively.  Thus,  we  may  readily 
determine  the  linear  combinations  of  the  components  of  tensors 
x„x,Xy,x,XyX|^, . . .  ,Sjj,SijS^„ . . .  ,5,jX*,5,jXnX,, . . .  which  belong  to  the  various 
irreducible  representations  yo,y,,y2,  •  •  • ,  F,,F2, . . .  of  the  group  F,. 

Let  us  consider  the  problem  of  determining  the  linear  stress-strain 
relation  for  a  material  whose  symmetry  is  defined  by  the  group  T,.  We  write 
the  constitutive  expression  as 

T  =  C,E,  (20) 

where 


T—  lit,  ,  +f22’ ^33’ ^1  l  ^22  +  ^‘^1  2’ ^1  1  ^22  2if,2, 


^13+'^23’^I3  '^2311^ 
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E,  =  (kn +^22.4^33.^11-^22  +  21^12,  Cl  1-^22-21^12, 

«13  +  '«23.«13-ie23r  (21) 

and  where  Ci  is  a  6  x  6  matrix.  If  we  refer  the  expression  to  the  reference 
frame  x,  whose  base  vectors  are  given  by  e,  =  Qij{6)ej,  we  have 

T  =  CiE„  T=R(0)T,  E,=R(0)E„  Ci  =  R(0)C,R-‘(0)  (22) 

If  the  reference  frame  x  is  an  equivalent  reference  frame,  we  require  that 
Cl  =Ci,  i.e. 

R(0)C,=C,R(0)  (23) 


We  observe  from  (14)  and  (18)  that 

R(0)=diag  (1,  1,  e~^‘*,  e^‘*,  e"‘*,  e*®)  (24) 


With  (23)  and  (24),  we  have  36  equations  relating  the  entries  C,j(j,y=  1, ...  6) 
of  Cl  which  are  given  by 


Cii=Cii,  Ci2  =  C,2,  Ci3  =  e-^‘®Ci.„  C.,  =  e^‘®C 

Ci5=e  ‘®Ci5,  Cig=e‘®Cig, . . . 

With  (25)  we  have 


Cii. 

Cl  2. 

0, 

0, 

0, 

0 

C21. 

C22. 

0, 

0, 

0, 

0 

0, 

0, 

C33. 

0, 

0, 

0 

0, 

0, 

0. 

C44, 

0, 

0 

0, 

0, 

0, 

0. 

C55. 

0 

0, 

0, 

0, 

0, 

0, 

^66 

14. 


(25) 


(26) 


This  tells  us  each  entry  in  T  which  belongs  to  a  representation  is 
expressible  as  a  linear  combination  of  the  elements  of  Ei  which  belong 
to  Vp.  Thus,  with  (21)  and  (26),  we  see  that  T  =  CiEi  may  be  written  as 


hi-^hz 

Cl  1,  C12 

^11+^22 

^33 

C21'  C22 

^33 

^11—^22  +2iti2  =C33(«’,  1—^22  +2i^l2).  7 2 
^11  +*22  — 2ill2  =  C44,(ei,  —^22  — 21^12).  ^2 
*13+i*23  =  C55(ei3  +  ie23)'  Vl 


*13  1*23-^66(^13  i4’23'('  ^1 
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where  the  Vo-  •  ■  ■  indicates  that  the  quantities  in  the  preceding  equation 

belong  to  the  irreducible  representation  yo . In  (27)  it  is  clear  that  we 

should  set  ^^=^33  and  If  we  set  C^^  =  a  +  ib,  C55  =  c  +  i<l,  the 

expressions  (27)3, . . .  ,(27)6  written  as 


Ill  I22 

a,-b 

^11  ^22 

Il3 

c,-d 

^13 

2tj2 

b,  a 

2e,2 

» 

123 

d,  c 

^23 

Let  us  consider  the  case  where  the  constitutive  expression  is  given  by 

^ij~^ijklmH^kl^mn’  ^kl~^lk  (29) 

We  write  this  in  matrix  form  as 


(30) 

where  T  is  given  by  (21)  and  Ej  denotes  the  (21  x  1)  column  matrix  whose 
entries  are  linearly  independent  linear  combinations  ot  ic  21  quantities 
el^,e^^e^2t■  ■  ■  so  chosen  that  each  belongs  to  one  of  the  irreducible 
representations  of  T^.  With  the  notation 

^2~^3J>  ^3  =  ^13'i''^23»  ^4"=^13~'^23 
£5  =  ^11 -«22  +  2i^l2.  £6  =  ^U-^22-2iC,2  ^ 

we  find  that  the  21  quantities  of  degree  2  in  the  Ei  which  belong  to 
yo,yi,ri,..,  are  given  by 

yo-Ei,  EiE2,  Eiy  E^E^, 


y,:E,E2,E2E2,E^E,; 
y2.E^E^,  £3^5’  ^3’ 
y^'EjEs'y 

'/V-Ei; 


r2.E,E^,E2E(„Ei,  (32) 

r^’-E^E^', 

r4-.£| 


The  constitutive  expression  T=C2E2  may  then  be  written  as 


I11  +  I22 

= 

^1’  *^2’  ^3>  ^4>  ^5 

^33 

Cft,  C-j,  Cg,  Cg,  CjQ 

£2 

E5E, 
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^13  '*■*^23  +<^12^2^3  +^13^4^5’  Vl 

tl3  —  it23=C,,E,E^+Ci2£2^‘i'^^I3^3^6’  ^1 
f  1 1  —  <22  "*■  2i<i  2  ~  ^14^1^5  +  ^15  ^2  ^5  ^1  6  ^3"  '/2 

^1  1  ~  ^22  “2*^12  ~^14^1^6  ^1  5^2^6  ■!■  ^16^41  Tj  (33) 
where  we  have  noted  that  £4  =  £3  and  E^  =  Ey 


3.  THE  [RREDUCIBLE  REPRESENTATIONS  FOR  THE 
TRANSVERSELY  ISOTROPIC  GROUPS 


There  are  five  groups  which  we  refer  to  as  transversely  isotropic  and  which 
are  denoted  by  T,.. . Tj.  These  groups  are  defined  by  listing  matrices  such 
that  these  matrices  or  products  of  these  matrices  specify  all  of  the  symmetry 
operations  associated  with  the  material  under  consideration.  We  may  refer 
to  these  matrices  as  generators  of  the  group.  We  then  define  the  irreducible 
representations  associated  with  a  group  7]  by  listing  the  matrices  which 
correspond  to  the  generators  of  the  group. 

Suppose  that  we  are  given  that  the  quantities  u,  .02.1103,041!  ^jlOj.Oft  1^ — 
belong  to  the  irreducible representaIions7o.‘;,,';2-73.-  •  ofagroupand  that 

the  quantities  h,.h2’il^3’KII^-  _  belong  to  the  irreducible 

representations  7o,7,, 72.73 _ of  the  same  group.  We  need  to  determine  the 

linear  combinations  Ci  =  a,jtajh^  of  the  products  of  the  Uj  and  which 
belong  to  the  various  irreducible  representations  of  the  group.  This 

information  is  provided  below  for  the  groups  T, . Tj  in  tables  which  are 

referred  to  as  product  tables.  Xu  ef  a]}  have  indicated  how  these  tables  may 
be  employed  in  conjunction  with  a  computer  program  to  automatically 
generate  constitutive  expressions.  The  extension  of  the  results  in  Ref.  2  to 
include  the  transversely  isotropic  materials  requires  the  development  of 
a  number  of  computer  programs.  This  work  will  be  carried  out  subse¬ 
quently. 

We  further  list  in  tables  entitled  basic  quantities  the  linear  combinations 
of  the  components  of  polar  vectors  p-.  axial  vectors  u,  and  symmetric 
second-order  tensors  which  belong  to  the  various  irreducible  re¬ 
presentations  of  the  group  considered. 


80 


G.  F.  Smith  and  G.  Bao 


3.1.  The  Group  Tj 

The  group  is  comprised  of  the  matrices  Q(0)  defined  by 


Q(0)= 


cos  6, 

sin  6, 

0 

—  sin  9, 

cos  9, 

0 

,0^9<2n 

(34) 

0. 

0, 

1 

The  group  T,  defines  the  symmetry  of  a  material  which  possesses  rotational 
symmetry  about  the  Xj  axis.  The  irreducible  representations  associated 
with  the  group  T,  are  ail  one-dimensional  and  are  given*  by 

yo-  1 

e-‘'^(p=l,2,....)  (35) 

r^:  e‘'’»(p=l,2....) 

In  (35)  the  1  X  1  matrices  1,  e"**^  and  e’'’*  correspond  to  the  group  element 
Q(0).  The  product  table  is  given  in  Table  1,  with  the  basic  quantities  in  Table 
2. 

3.2.  The  Group  Tj 

The  group  T2  is  comprised  of  the  matrices  Q{6)  and  R,Q(0)  where 


Table  1 

Product  table  for  7, 


Vo-  ^0’  ^0 

<^0^0 

. ) 

yp' 

tto^p’  ^O^p 

=  1,2, . . . . ;  m  +  n  =  p) 

OmB,,  A„b„  (m,n=l,2,...;  m-n  =  p) 

r.:  A„ 

P  P*  P 

t^oBfp  ApbQ 

A„B„  (w,n  =  l,2,...;  m  +  n  =  p) 

A„b„,  a,B„  (m,n  =  1 ,2, . . . ;  m  -  n  =  p) 
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Table  2 

Basic  quantities  for  T, 


Vo-  P3>  ^22>  ^33 

vr-  Pi+ip2,ai+iaj,S,3  +  iS23 

rp  Pi  ~*P2>  ~'^2' ^13  “‘^23 

yi-  ^11  ~  ^22  ■*■2*^1 2 

r2’  ^U~^22~2iSj3 


O^0<2;t  and  where 


W)= 


cosfl. 

sin  0y 

0 

—  sin  0,- 

cosd. 

0 

,R,=diag(- 1,1,1) 

(36) 

0, 

0, 

1 

The  irreducible  representations  associated  with  the  group  Tj  are  defined  by 
listing  the  matrices  corresponding  to  the  group  elements  Q(0)  and  R,.  We 
denote  the  irreducible  representations  (see  Ref.  5)  by 


ro:  >.  1 


ro-.  1,-1 


e'*'^,  0 

l|0,  1| 

0,  e**^ 

O 

(37) 


Table  3 

Product  table  for  Tj 


Vo-  ^0' 

'^0^0 

0«l^’».2+“».2<’ml  (W=l,2,....) 

r  o'  '^0’ 

'^0^0 

(m=l,2 . ) 

iiapi,s2iiMi^i,^2ir 

^0^p2ll^,  ll^plf’o,  ^(72boll 
IMof’pl,  ~  ll<*plSo,  — ^p2®oll^ 

^mibp2  r.  (m,  n  =  1, 2 . ;m  +  n  =  p), 

ll‘»mi^2,  a».2<’,i  r.  Ila,2<’».l,«l.l*«.2ll^ ('«,"=  1. 2 . ;  m-n  =  p) 


V 
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Table  4 

Basic  quantities  for 


Vo: 

Pi’  '^^22’  ^33 

W 

O3 

Vv 

llPi-hiPi,  -p,+ip2ll^’  l|a,+W2. 

11^13  ~  3  "*"*^2311 

Ji- 

|lS,,-S22  +  2iS,2,S,,-S3j-2iS,2r 

where  the  first  and  second  matrices  correspond  to  Q(0)  and  R,  respectively. 
The  product  table  is  given  in  Table  3  and  basic  quantities  in  Table  4. 


3.3.  The  Group 

The  group  T3  is  comprised  of  the  matrices  Q(0)  and  R3Q(0)  where 
0^6 <2n  and  where 

cos  6,  sin  6,  0 

Q(0)=  -sin0,  cose,  0  ,  R3  =  diag(l,l,- 1)  (38) 

0,  0,  1 

The  irreducible  representations  associated  with  the  group  Tj  are  defined  by 
listing  the  matrices  corresponding  to  the  group  elements  Q(0)  and  Rj.  We 
denote  the  irreducible  representations  by 

To:  1,-1 

V  e--'’M;ve-'^,l  (p=l,2,..,.)  (39) 

r^:  e-‘'^,-l;rp:e-'^, -1  (p=1.2 . ) 

where  the  first  and  second  1  x  1  matrices  correspond  to  Q(0)  and  R3 
respectively.  The  product  table  is  given  in  Table  5  and  basic  quantities  in 
Table  6. 

3.4.  The  Group 

The  group  T4  is  comprised  of  the  matrices  Q(0),  R,Q(0),  RsQf^)  and 
R,R3Q(0)  where  0^6<2ti  and  where 

cos  6,  sin  6,  0 

Q(0)=  -sin0,  COS0  0  ,  R,  =diag(- 1,1,1),  R3  =  diag(l, 1,-1)  (40) 

0,  0,  1 
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Table  5 

Product  table,  Tj 


^0'  ^0 
Aq^O'  ^0^0 

(m=  1,  2 . ) 

■^0’  ®0 
^0®0>  ^^0^0 

amS„,  a„B„,  (m  =  1, 2, . . . . ) 

ttobp,  apbQ,  AqBp,  ApB^ 

(m,n  =  \,2,...;m  +  n  =  p) 

^A.  (m,  n=l.2,...:  m-n  =  p) 

dp,  ^p 

flfl^p’  ^p^O’  '^0®p’ 

dm^,v  (m,n=l,2,^..:m  +  n  =  p) 

dj,,  <jA.  ABp.  -4,3^  (m,  n=l,  2 . ;  m-n  =  p) 

-4p,Bp; 

tt(yBpy  v4pbQ,  A^bp,  QpBff 

A„Bp{m,n=  1,2 . ;  m  +  n  =  p) 

a„S„,  d,B„,  (m,  n  =  1 , 2 . ;  m  -  n  =  p) 

/l„g„ 

P  P 

aoBp,Apbo,AoBp,dpB(, 

dmB„,  -4A„,(m,n  =  l,2,....;m  +  n  =  p) 

0™®,.  a,g»P  (m,  n  =  1. 2, . . . ;  m  -  n  =  p) 


Table  6 

Basic  quantities,  Tj 


Vo-  ‘•3’^11+'S22'533 

^0'  P} 

Pl+iP2 
Vi:  Pi->P2 

Fi.'  a, +ia2, 5(3 +iS23 

Fp  a,  —  ia2,  5|3  —  iS23 

yji  Sii  — S22  +  2iS,2 

V2'  ^1  1  “^22  ~  2iS,2 
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The  irreducible  representations  associated  with  the  group  are  defined  by 
listing  the  matrices  corresponding  to  the  group  elements  Q(0),  Rj  and  Rj. 
We  denote  the  irreducible  representations  by 


Vor 

1,  1, 

1 

yo2- 

1,  1,- 

1 

703-- 

1,-1, 

1 

yo4- 

1,-1,- 

1 

(41) 

e-ip«^ 

0 

0, 

1 

1,  0 

yp- 

♦ 

(p=l,2,....) 

0, 

gip« 

1, 

0 

0,  1 

e-ip«^ 

0 

0. 

1 

-1, 

0 

rp- 

(p=l,2,...) 

0, 

gipe 

1, 

0 

0, 

-1 

where  the  first,  second  and  third  matrices  correspond  to  Q(0),  R,  and  Rj 
respectively.  The  product  table  is  shown  in  Table  7  and  basic  quantities  in 
Table  8. 


3.5.  The  Group  Tj 

The  group  Tj  is  comprised  of  the  matrices  Q(0)  and  DjOfS)  where 
0^0 <2n  and  where 

cos  0,  sin  0,  0 

Q(0)=  —  sin0,  COS0,  0  i,D2  =  diag(-l,  1,  -1)  (42) 

0,  0,  1 

The  irreducible  representations  associated  with  the  group  Tj  are  defined  by 
listing  the  matrices  corresponding  to  the  group  elements  Q(0)  and  Dj.  We 
denote  the  irreducible  representations  by 

n.:  I, -1 


o:  I.-l 

0  0,  1 
'p- 

n  »'p»  1  n 


(/>=!, 2,....) 
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Table  7 
Product  table  for 

voi;  oi.i’i 

d^b ^2^2* 

(m=l,2,...) 

y02‘  <^2^^2 

^\b2y  ^2^1a^3^4.a^4.^3 

^ml^m2  ^m2^ml*  ^mi^n2  ■^m2^ml  (m=1.2,...) 

Vos- 

a^b^,a2b^,a^b^,a^b2 

"^ml  ®m2  ~ (W=  1,  2,  .  .  .  .  ) 

V04:  ‘»4.  K 

^1^4»  ^2^3*  ^3^2*  ^4^1 

‘^ml®i«2  ~^m2®iiil>  '^ml^iii2  ~ '^iii2^ml  {W=  1,  2,  .  .  .  .  ) 

Vr  I|a,i,ap2ll^ll*pi.^p2r 

ii^i^pf‘*i*p2ir,  ii«pi<’i.‘*p2*ir 
ll^2®pl'  ^2®p2ll^'  Mpl^2’  ^p2^2ll^ 
ll‘J3*pl.  -a3^p2r.  Il“pl*»3.  -‘>p2^’3ll’' 

Mp,fc4.  -^p2Mr 


II  1  1  >  ^m2  ^«2  II  II  -^ntl  Bp  1  >  '4|n2  ®»2  II  ^ 

(m,  n  =  1 , 2, . . 

.  ;  227  +  27  =p) 

l|amlfr,2.  «»>2fr,l  ir,  ll«»2i’ml.  Opl*>m2ir 

(m, «  =  1 , 2, _ 

;m-n  =  p) 

M„lB,2M„2Bpir,  M,2Sp..M„.B„2ll 

^  (m,n=  1,2,, 

...■,m-n  =  p) 

Fp:  Mp„/lp2ir,IIBp„Bp2ir 

ll^i^pi’  ‘*i®p2ll^>  Mpi^n  ‘4p2f’i  11*^ 
ll^2^pl>  ^2^p2ll  ll^pl^2’  ®p2^2ll^ 

ll^3®pl>  ~^3®p2ll^>  ll'4piB3,  “-41,2^311^ 
ll"4^’pl.  -“4<’p2r.  Il«pl<>4.  -"p2^’4r 

ll^iiil®nl>  ‘*ifi2®»2ll^'  ll'4»,)l>,i.  '4„2i’p2ll^ 

(m,  21=  1,2,.. 

.;2n  +  27  =  p) 

ll“iiil®«2>  ‘*iii2®»I  II  Il<*p2®iiil’‘*»l®m2ll^ 

(m,  27=  1, 2, . . 

.;m-n  =  p) 

m„.6.2. '4p.2<’p.  ir.  Mp2^..  ^..<’«2ir 

(227,  27=  1,  2,  .  . 

.;  277  -  27  =  p) 

where  the  first  and  second  matrices  correspond  to  Q(0)  and  Dj  respectively. 
We  note  that  the  irreducible  representations  (43)  are  the  same  as  those 
appearing  in  (37).  The  product  table  will  then  be  the  same  as  the  product 
table  for  which  is  shown  in  Table  9,  with  basic  quantities  in  Table  10. 
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Table  8 

Basic  quantities  for 

yoi- 

^1  1  +^22'  ^33 

7o2- 

Pi 

>’03- 

tti 

)'04- 

Vl 

IIPi+iP2.  -Pi+iPiW^ 

r.: 

llai+iaj.ai-ia^r,  IlSu  +  iSjj,  -Su  +  iS^jir 

Va- 

||S,,-S22  +  2iS,2,S,.-S22-2iS,2r 

Table  9 

Product  table  for  T, 

Tp-  ^0'  ^0 

^0^0>  '^0^0 

r  o'.  Aq,  Bg 

^0®0'  '^0^0 

(m=l,2 . ) 

V  ii«p,.«p2iiMi^..^2ir 

-/lobp2ir.llapiBo. -«p2Bor 

II  "o^p  1 .  «0^p2  II  Ml  "p  1  ^0-  «p2^0  II  ^ 

I|ami^i.am2bp2ir,  («,  21=1,2 . ;m  +  n  =  p) 

Il0p,i^2. OmiK ,  II ^  lla,2«’».i’ "-I^m2 (m,  H  =  1 , 2 . ;m-n  =  p) 


Table  10 
Basic  quantities.  T, 

Vo-  1 +^22>  ^33 

Tq-  Ps,  <*3 

>’i:  IIPl+'P2.  -Pi+>P2ir.  ll«l+'«2.  -«l+ia2r-  l|S,3  +  iS23.'Sl3-'S23ll^ 

yi'  ll^i  I  “^22  2iS|2,  S, ,  —  Sjj  —  2iS,2l| 
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4.  APPLICATIONS 

In  this  section,  we  give  some  examples  of  the  application  of  the  results 
derived  above  to  the  generation  of  non-linear  constitutive  expressions.  We 
first  consider  the  problem  of  determining  the  form  of  a  second-order 
tensor-valued  function 


Tij=Cij^t„x^x,x„  (44) 

which  is  of  degree  three  in  the  components  of  a  polar  vector  Xj  and  which  is 
invariant  under  the  group  T,.  From  Table  2  we  see  that 

^1 1  “b  ^22>  ^33»  ^13  “b  *^23’  ^13  “*^23'  ^1 1  ~  ^22  "b  2it  j2,  t ]  1  —  t22  ~  ^it j2  (45) 

belong  to  y^,  ’/q,  y,,  Tj,  y2,  r2  respectively  and  that 

X3,x, -|-ix2,  X,  —  ix2  (46) 

belong  to  yo.yi.ri  respectively.  Upon  employing  Table  1  twice,  we  see  that 

7cix3(x?  +  xi),x|(xi-l-ix2),(xi  +xi)(Xi  4- ix2),xi(Xi -ix2),(x? +  xi)(x, -^2), 

X3(xf-X2  +  2ixiX2),  Xjlxi -xf  — 2ix,X2),(x,  -l-ix2)^,  (Xi  -ix2)^  (47) 

belong  to  yo,yo.yi.yi>ri.ri.y2T2.73.r3  respectively.  Each  of  the 
quantities  in  (45)  which  belongs  to  a  representation  y^  (say)  is  expressible  as 
a  linear  combination  of  the  quantities  in  (46)  which  belong  to  -/p.  Thus,  we 
have 


hi+hz 

Cp  C2 

hi 

^■3'  ^4 

X3(X?-|-xi) 

^  3  +  '^23  =  <^5^3(Xi  +  1X2)  +  C^(xf  +  xfXXj  +  IX^} 

fl3  -‘f23=^5^3(7Ci  -iX2)-|-C6(x?  -I-X^XX,  -iX2)  (48) 

fll-^22  +  2iti2=‘‘77C3(^f -7^2  +  21X1X2) 
tu-f22-2itl2=M3(7Cl-7C2-2ix,X2) 

where  C5, . . . ,  c,  denote  the  complex  conjugates  of  C5, . . . ,  C7. 

We  next  consider  the  problem  of  determining  the  form  of  the  polar 
vector-valued  function 


X-,  —  +  ^ijklm^jk^tm 


(49) 


of  the  symmetric  second-order  tensor  S,y  which  is  invariant  under  the  group 
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Tj.  We  employ  the  notation 

■^1  ~  “^1 1  ^2  ~^33»  ^3  ~^I3  "i"  ’^23'  ^4  ~  ~  ^13  "i"  ’^23 

^5  ”■^22  +  2iS|2)  ^6~^ll~^22 

We  note  that  S^  =  —S^  and  —  S^.  From  Table  4  we  see  that 

X3,  ||x,+ix2, -Xj+ixjl!^  (51) 


belong  to  and  Vj  respectively  and  that 


s„s2,  iis3,s4ir.iis5,sgir 

(52) 

belong  to  Vq)  yoi  ^nd  y2  respectively.  In  (52)  we  have  used  the  notation 
(50).  Upon  employing  Table  3,  we  see  that  the  21  products  of  degree  two  in 
the  Si(i=  1, _ ,6)  belong  to  the  representations  listed  below. 

To 

5i,  si,  SjS2,  S3S4,  SjSg 

?! 

s,  IIS3,  S4II ^  S2IIS3,  S4II  ^  IIS4S3,  S3S,  II  ^ 

72 

s,iis5.,s,r,s2iis5.s6irjisisir 

(53) 

73 

l|S3S5,S4Sjr 

74 

((si.sir 

The  constitutive  equation  (49)  may  then  be  written  as 

X3=CiSi  +C2S2  +  C3S1 +C4S2  +  C5S,S2 +  0^X354 +  C7S5S6 


X,  -t-iX2 

_ 

Cg,  0 

S3 

-1- 

C9, 0 

S.S3 

+ 

0 

0 

S2S3 

-X,  -|-iX2 

O.Cg 

S4 

0,  Cg 

S,S4 

1  0,  c,o 

S2S4 

Cij,  0 

S4S5 

0,  Cji 

S3Se 

(54) 


Since  — x, +  ix2  =  — (x, +  ix2)  and  S^=-Sj,  we  see  that  c^  =  c^  which 
implies  that  Cg  is  a  real  number.  Similarly,  we  see  that  Cg,  Ck,  and  Cj  1  are  also 
real  numbers. 

We  next  consider  the  problem  of  determining  the  form  of  the  polar 
vector-valued  function 


~ ^ijk^Jk  ^ijklm^jk^lm  (55) 

and  the  axial  vector-valued  function 


~  ^ijk^jk  ‘^ijklm^jk^lm 


(56) 
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of  the  symmetric  second-order  tensor  S^j  which  is  invariant  under  the  group 
Ty  We  employ  the  notation 


1  ^2  ^3  ~  ^13  ■i'*^23>  ^4  ~  ^13  *^23 

~^22  +  2i5j2.  ~'^22  ~2i5i2 

We  note  that  S^  =  Si  and  S(,  =  Sy  From  Table  6  we  see  that 

X3,  X,  -f-ix2,  Xi  —  ix2 

belong  to  Fq,  yi.y,  respectively,  that 

flj.fl,  +ia2. «!  —  i^i2 

belong  to  J'l*  respectively  and  that 


(57) 


(58) 

(59) 


^It  ^2.  Sj,  S5,  Sg 


(60) 


belong  to  yo.  ?o>  I'l*  ^i’  >’2>  72  respectively.  Upon  employing  Table  5,  we  see 
that  the  21  products  of  degree  2  in  the  S((i=  1, _ ,6)  belong  to  the  re¬ 

presentations  listed  below. 

Vo-  ^2>  ^1^2>  ^3^4>  ^5^6 


F,:  S1S3,  S2S3,  S4S5 
Fj;  SjS4,  S2S4,  SjSg 
72-  ^1^5’  ^2^5’  ^3 
V2:  S,S„S2S„Si 

^3-  SySs 
^3-  S4S. 

74:  SI 
74:  Si 


(61) 


With  (58), . . .  (61),  we  see  that  there  are  no  terms  of  degrees  1  or  2  in  the 
Sj(i=  1, . . .  ,6)  which  belong  to  any  of  the  representations  to  which 
X3,  Xi  +ix2,  X,  —  ix2  belong.  This  implies  that  a  constitutive  expression  of 
the  form  (55)  is  ruled  out  by  symmetry  considerations.  The  constitutive 
expression  (56)  may  be  written  as 

~  d iS I  -t- d2S2  "i- d^Si  -i-  d^Sl  -h  d^S  1S2  d^S^S^  d-jS 
a^  -h  102  ~d^S2  +  dgS^S2  +  <^io'52'53  +d^  iS^S^ 

Uj  —ia2  =  SgS^  +  3gS^S^  +  3^QS2S^  +  3^^S2S^ 


(62) 
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With  the  notation  <i8=^8  +  */8» _ *  +i/ii>  we  see  from  (57)  that 

(62)2  3  ^  written  as 


®8>  “/s 

^lal 

1 

-1-  • 

■  -1- 

/ll 

(Sn 

~  ^22)^1 3 +  2Sj  2^23 

^2 

/s’  ^8 

^23! 

/u»  ^11 

-(Stt 

~  522)^23 +  2S12S1 3  1 

(63) 
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ABSTRACT 

This  paper  focusses  on  the  study  of transverse  cracking  in  cross-ply  composite 
laminates  under  monotonic  load.  In  particular  an  explicit  formula  is  obtained 
for  the  loss  of  stiffness  of  a  cracked  laminate.  But  the  major  part  of  the  paper  is 
devoted  to  the  statistical  fracture  mechanics  of  progressive  cracking.  We  show 
how  to  predict  crack  density  as  a  function  of  applied  load.  The  model  agrees 
well  with  experiment. 


1.  INTRODUCTION 

Problems  associated  with  transverse  matrix  cracking  in  cross-ply  com¬ 
posite  laminates  have  been  discussed  by  many  investigators.  Two  partic¬ 
ular  goals  have  been  to  predict  loss  of  stiffness  and  crack  density  under 
monotonic  tensile  load.  To  the  best  of  our  knowledge  the  problem  was  first 
addressed  by  Garrett  and  Bailey.*  The  model  proposed  by  these  authors  is 
usually  referred  to  as  shear  lag  theory.  This  model  was  subject  to  further 
refinement  and  development,  culminating  in  the  paper  by  Bailey  et  al.^ 
Further  contributions  to  the  study  of  transverse  cracking  have  been  given 
by  a  number  of  authors.  For  a  reasonably  complete  account  see  Laws  and 
Dvorak.^ 
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The  aim  of  the  present  paper  is  to  give  a  succinct  account  of  the 
progressive  transverse  cracking  theory  developed  elsewhere  by  the  present 
authors.  The  mechanical  model  described  is,  perhaps,  the  simplest  one¬ 
dimensional  model  possible.  Indeed,  some  straightforward  stress  analysis 
enables  us  to  predict  the  loss  of  stiffness  for  a  given  crack  density.  This 
prediction  compares  well  with  experiment^  and  with  the  self-consistent 
model  of  Laws  and  Dvorak"*  and  the  lower  bound  given  by  Hashin.^ 
However,  from  a  practical  point  of  view,  the  crucial  problem  is  to  predict 
the  transverse  crack  density  as  a  function  of  the  applied  (monotonic)  load. 
This  problem  is  solved  herein  by  first  analysing  the  fracture  mechanics  of 
transverse  cracking  and  second,  by  giving  an  analysis  of  the  statistics  of 
progressive  cracking.  The  result  is  a  formula  which  permits  us  to  predict 
transverse  crack  density  as  a  function  of  applied  load.  Theory  is  shown  to 
compare  well  with  experiment. 


2.  PRELIMINARIES 

In  this  paper  we  consider  monotonic  tensile  loading  of  cross-ply  composite 
laminates,  see  Fig.  1.  It  is  convenient  to  use  the  subscripts  t  and  1  to  refer  to 
the  transverse  and  longitudinal  plies  respectively.  Thus  the  initial  stresses  in 
the  laminate  are  respectively  of  and  of. 

>b 


Fig.  1.  Symmetric  cross-ply  laminate  under  axial  load. 


The  one-dimensional  theory  which  is  outlined  here  assumes  that  the 
displacement  of  each  layer  is  constant  over  the  thickness  of  that  layer.  As 
usual,  the  displacements  u(x)  and  v(x)  are  measured  from  the  state  in 
which  there  is  no  mechanical  loading.  The  associated  strains  are 

dp  du 


0) 
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When  the  laminate  is  subject  to  mechanical  loads,  the  total  stress  is  the  sum 
of  the  residual  stress  plus  the  stress  due  to  loading: 

ff,  =  <Tf  +  T„  T,  =  £,e,  (2) 

ff,  =  fff  +  T„  T,  =  £,£,  (3) 

where  E  denotes  Young’s  modulus. 

Turning  now  to  the  equilibrium  of  the  laminate  it  is  clear  that 

baf  +  daf=0  (4) 

and  that  under  applied  stress 

ba,+da,=(b  +  d)(7^  (5) 

The  essential  feature  of  shear  lag  theory  is  that  the  shear  stress  which  is 
responsible  for  the  load  transfer  between  the  0°  and  90°  plies  is  proportional 
to  the  relative  displacement  of  the  two  layers.  In  other  words 

z  =  K(v-u)  (6) 


where  K  is  a  constant. 

Finally  we  recall  from  Ref  3  that  equilibrium  of  the  0°  plies  requires  that 


T  = 


(7) 


whereas  equilibrium  of  the  90°  plies  demands  that 


dx 


(8) 


Within  the  framework  of  the  theory  developed  here,  Young’s  modulus  for 
the  uncracked  laminate  is  given  by 


bE,+dE^ 


(9) 


Stress  analysis  of  the  cracked  laminate  is  most  easily  achieved  with  the 
help  of  the  differential  equation  for  the  transverse  stress.  This  equation  is 
obtained  by  differentiating  eqn.  (8)  and  making  use  of  (l)-(6)  together  with 
(9)  to  give 


Here,  we  have  introduced  the  non-dimensional  shear  lag  parameter 
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^  through 

Kd(bE,+dE,) 

^  b£,£,  '  ’ 

We  now  focus  our  attention  on  the  ligament  between  adjacent  transverse 
cracks,  see  Fig.  2.  In  this  configuration,  the  differential  equation  (10)  holds 
in  the  ligament  AB  and  the  correct  boundary  conditions  for  eqn.  (10)  are 

a^=0  when  x=±h  (12) 


- 2d/(3 - ► 

Fig.  2.  Two  adjacent  transverse  cracks  in  the  90‘  plies. 


The  solution  for  the  transverse  stress  is 


As  is  shown  by  Laws  and  Dvorak,^  it  is  not  difficult  to  evaluate  the 
displacements  of  the  respective  plies.  For  our  present  purposes  it  suffices  to 
note  that  the  average  strain,  £„  of  the  ligament  AB,  as  measured  at  the 
surface  of  the  laminate  is 


d^E, 

^hbE, 


+ 


d^(7f 

^hbE, 


(14) 


We  see  from  (14)  that  the  ligament  has  acquired  a  permanent  strain, 
given  by 
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Clearly  the  permanent  strain  (15)  is  due  to  initial  stress.  However,  it  turns 
out  that  for  all  practical  problems  of  interest,  the  value  of  e^,  is  insignificant 
compared  with  e,.  This  simple  theoretical  argument  supports  the  observa¬ 
tion  that  any  permanent  strain,  due  to  the  relaxation  of  initial  stress  at  crack 
surfaces,  is  negligible.  Now,  for  a  cracked  laminate,  in  which  the  transverse 
crack  density  is  the  average  distance  between  cracks  is,  from  Fig.  2, 

2h  =  2d/P 

Hence  from  (14)  we  obtain  the  formula  for  the  Young’s  modulus  £(^),  of  the 
cracked  laminate: 


£(«-£.{l  +  |^lanh5|  il6) 

It  is  clear  from  (16)  that  in  order  to  predict  the  loss  of  stiffness  of 
a  cross-ply  laminate  we  need  to  know  the  value  of  the  shear  lag  parameter 
This  has  been  argued  in  depth  by  Laws  and  Dvorak.^  Briefly,  the  suggestion 
is  that  ^  be  determined  by  measurement  of  the  first  ply  failure  stress  —  as  is 
discussed  shortly.  It  turns  out  that  the  loss  of  stiffness  is  not  particularly 
sensitive  to  the  value  of  It  is  also  appropriate  to  mention  that  the 
self-consistent  model  of  Laws  and  Dvorak,"*  the  lower  bound  obtained  by 
Hashin,*  as  well  as  the  shear  lag  modeP  give  excellent  agreement  with 
experiment.  For  details  we  refer  the  interested  reader  to  Ref.  3. 

In  this  contribution  we  prefer  to  focus  on  the  more  difficult  problem  of 
predicting  the  transverse  crack  density  (fi)  as  a  function  of  the  applied  load 
(<r,).  Once  this  relationship  is  available  it  is  clear  that  we  can  also  predict  the 
appropriate  stress-strain  relation,  etc. 


3.  PROGRESSIVE  CRACKING 

The  detailed  analysis  of  progressive  transverse  cracking  is  given  by  Laws 
and  Dvorak^  so  we  are  content  here  to  emphasize  the  main  results. 
Consider  the  uncracked  ligament  of  Fig.  2.  When  the  applied  load  reaches 
a  critical  value  this  ligament  will  itself  crack  at  some  location  C,  as  in  Fig.  3. 
Since  the  transverse  cracking  process  is  not  deterministic  there  is  no  reason 
to  suppose  that  C  lies  at  the  mid-point  of  AB.  It  therefore  follows  that  we 
need  to  compute  the  energy  release  rate  for  a  crack  at  an  arbitrary  location 
C  propagating  across  the  ply.  The  required  calculation  is  not  easy. 
Nevertheless  we  can  read  off  the  required  result  from  the  paper  of  Laws  and 
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A  C _ ^ 


hg - ► 

- 2h - ► 

Fig.  3.  The  ligament  AB  with  an  additional  crack  at  C. 


Dvorak:^ 


d{b  +  d^Eo 


tanh  ^  +  tanh  ^  -  tanh 
2d  2d 


?} 


(17) 


To  get  the  required  energy  release  rate  at  first  ply  failure  we  consider  the 
limit  in  which  h,  It,  and  /ij-^oo,  namely 


Qtvt 

^bE,E, 


Clearly  we  get  first  ply  failure  when 


Thus  it  follows  that 


d(b+d)Eo 

G,bE,E, 


(18) 


Now  standard  data  for  cross-ply  laminates  is  given  by  Wang®  for  two 
different  graphite-epoxy  systems.  In  short,  all  quantities  in  eqn.  (18),  except 
(J,  are  either  given  or  easily  inferred.  Thus  we  regard  eqn.  (18)  as  the  rule 
which  determines  the  shear  lag  parameter  Values  for  specific  systems  are 
given  later. 

We  now  turn  to  the  problem  of  determining  crack  density  under 
increasing  load.  Here  it  is  necessary  to  make  some  assumptions  about  the 
statistics  of  progressive  cracking.  In  this  connection  we  refer  to  some 
arguments  of  Laws  and  Dvorak^  who  appealed  to  elementary  fracture 
mechanics  to  show  that  an  appropriate  choice  for  the  probability  density 
function,  p,  for  the  site  of  the  next  crack  in  the  ligament  of  Fig.  3,  is  that  p  be 
proportional  to  the  stress  in  the  transverse  ply.  Of  course,  normalization 
gives  the  factor  of  proportionality.  Thus  from  eqn.  (13)  with  x  replaced  by 
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(y  — #i)  we  find 


p(y)= 


cosh  — 
a 


tanh  — 
a 


It  now  follows  that  in  a  laminate  which  already  contains  cracks  of  den¬ 
sity  the  expected  value  of  the  applied  stress  to  cause  additional  cracking  is 


£K(/3))= 


p(y)<7,(y)dy 


The  integral  in  eqn.  (20)  must  be  evaluated  numerically  and  demands 
considerable  care. 

We  now  compare  the  theoretical  predictions  with  some  experimental 
results  of  Crossman  et  al.^  for  two  graphite-epoxy  systems.  Perhaps  the 
most  convenient  sources  for  the  data  are  the  survey  article  by  Wang® 
together  with  the  report.®  Use  of  these  data  enables  us  to  calculate  the 
appropriate  value  of  (  directly  from  eqn.  (18).  The  results  are,  cf.  Ref.  3, 


90), 

^  =  0-93 

90^), 

1) 

00 

*2> 

9O3). 

^  =  2-24 

Comparisons  of  the  theoretical  predictions  with  experimental  results  are 
shown  in  Fig.  4.  For  clarity  we  have  not  included  the  predictions  of  the 
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Applied  load  (MPa) 

Fig.  4.  Theory  versus  experiment  for  progressive  cracking  of  AS-3501-06  laminates. 
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Wang-Crossman  theory  in  Fig.  4.  However,  we  note  that  the  theory  of  these 
authors  also  gives  excellent  agreement. 

Finally,  we  turn  to  some  experimental  data  obtained  by  Wang®  for  some 
T300/934  laminates.  For  these  laminates  we  infer  the  values  of  to  be  as 
follows:^ 


(0, 

90^, 

0) 

=  108 

(0, 

9O3, 

0) 

=  1-70 

(0, 

90,^, 

0) 

=  1-79 

Theory  is  compared  with  experiment  in  Fig.  5.  Agreement  is  good. 


Fig  5  Theory  versus  experiment  for  progressive  cracking  of  7300/934  laminates. 
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ABSTRACT 

Predicting  failure  modes  and  mechanisms  of  failure  in  composites  hinges  to 
a  great  extent  on  our  knowledge  of  the  stress  fields  associated  with 
discontinuities  in  composites.  In  composites,  singularities  arise  at  discontinui¬ 
ties  at  the  micro  and  macro  scale  of  their  behavior,  and  whence  affect  the 
failure  mode  and  damage  at  all  levels.  Such  discontinuities  include  interfaces, 
boundaries  of  laminated  composites,  cracks  between  layers,  design  discontinu¬ 
ities  such  as  joints,  cut  outs,  cracks,  repairs,  etc.  The  singularities  in  many  of 
these  cases  are  not  known,  or  the  known  solutions  are  too  complex,  and  no 
general  method  is  at  the  disposal  of  the  designer  or  experimentalist  for 
analysing  the  specific  case  he  is  dealing  with. 

Recently,  the  author  developed  a  general  Finite  Element  based  iterative 
method  for  the  solution  of  eigenvalue  problems  common  in  the  fracture  of 
composite  materials.  The  advantage  of  the  method  is  that  it  relies  on  the  use  of 
general  purpose  finite  element  packages  for  performing  the  iterative  analysis. 
The  data  are  processed  to  evaluate  both  power  singularities  and  oscillating 
singularities,  occurring  at  interface  cracks  of  composites.  In  this  presentation, 
we  will  illustrate  several  problems  of  interface  cracks,  and  delamination,  and 
show  how  the  iterative  finite  element  method  is  used  in  evaluating  the  singular 
field.  The  implications  of  the  singularities  to  failure  mechanisms  of  fiber/ 
matrix  debonding,  edge  delamination  and  the  behavior  of  laminated  compo¬ 
sites  will  be  discussed. 
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1.  SINGULARITY  FIELDS  IN  COMPOSITE  MATERIALS 

The  asymptotic  field  in  many  singularity  problems  associated  with  the 
fracture  of  composite  materials  can  be  written  as* 

(la) 

and 

F}(e,  C,)+  ^r*-'f  ?(e,  Q)  (lb) 

where  <5  is  the  power  of  the  singularity,  which  could  be  real  or  complex 
depending  on  the  nature  of  the  discontinuity  and  material  properties.  The 
functions /"y  and  F"  are  also  dependent  on  the  materials  properties  C*  and 
are  both  bounded  functions.  The  restrictions  on  the  power  of  singularity  are 
that  it  produces  finite  displacements  at  the  origin  and  that  the  resulting 
strain  energy  is  positive  definite. 

LI.  The  Finite  Element  Iterative  Method 

The  aim  of  the  finite  element  iterative  approach  discussed  here  is  the 
evaluation  of  singularity  5  and  the  functionsXy(0)  and  f  ,(0)  in  eqn.  ( 1 ).  It  will 
also  be  shown  that  the  approach  is  a  global  method  and  whence  the  stress 
intensities  can  be  calculated.  This  technique  will  be  referred  to  as  the  Finite 
Element  Iterative  Method  (FEIM).^”’  Since  the  method  uses  a  displace¬ 
ment  formulation  with  elements  that  satisfy  all  the  requirements  of  rigid 
body  motion,  constant  strain,  positive  definite  strain  energy,  and  finite 
displacement  for  linear  elastic  problems,  the  conditions  on  the  asymptotic 
eigenvalue  solution  are  satisfied  a  priori. 

l.i.I.  Procedure 

FEIM  uses  existing  general  purpose  Finite  Element  Programs  as  follows: 

1.  Construct  a  circular  mesh  around  the  crack  or  singularity,  Fig.  1. 

2.  The  radii  of  the  rings  should  follow  an  (r^)  distribution.  For  crack 
problems  the  quarter-point  elements  are  used  to  represent  the  square 
root  singularity. 

3.  Impose  arbitrary  boundary  displacements  on  the  outer  boundary 

— an  approximate  solution  will  accelerate  the  convergence. 

4.  Choose  radius  R,  near  the  crack  tip  or  singularity,  where  the 
displacements  {u„_}  are  evaluated. 
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5.  Perform  a  finite  element  analysis  and  extract  the  displacements 
Scale  the  results  at  R,  and  apply  these  as  boundary  displacements 

at  the  outer  boundary 

6.  Repeat  step  5  until  convergence  is  attained. 

7.  To  improve  the  results,  subtract  the  crack  tip  displacements  u„  from 
the  displacements  at  the  rind  R,  before  scaling.  The  displacements  u„ 
represent  a  rigid  body  motion  and  therefore  do  not  affect  the  solution. 

The  above  procedure  was  originally  applied  to  the  case  of  a  crack  in 
a  general  anisotropic  material^  and  was  found  to  converge  to  the 


104 


R.  S.  Barsoum 


asymptotic  field  exact  solution.®  The  singularity  in  this  case  is  y/r,  and  the 
only  unknowns  are  the  distribution  functions  and 

Convergence  is  judged  when  the  displacements  at  any  ring  eventually 
approach  ~*FXQ)-  Hence,  5  and  Fi(9)  can  be  extracted  from  the  analysis. 

In  the  cases  where  the  singularity  5  is  a  real  number,  we  find  that  the 
displacements  }  do  not  change  from  one  iteration  to  the  next,  and  the 
scaling  parameter  reaches  a  constant  value.  The  asymptotic  field  in  this  case 
is  called  self-similar.  In  the  case  where  the  singularity  3  is  a  complex 
member,  the  displacements  have  a  phase  shift  from  one  iteration  to 
the  next  and  the  scaling  parameter  could  oscillate  wildly  from  one  iteration 
to  the  next.  However,  that  is  the  nature  of  the  problem  and  the  evaluation  of 
this  field,  which  is  called  a  non-self-similar  field,  will  be  discussed  later. 


2.  IDENTIFICATION  OF  THE  FEIM  WITH  THE 
POWER  METHOD  FOR  EIGENVALUE  PROBLEMS 


The  Finite  Element  Iterative  Method  was  identified  with  the  Power  Sweep 
method  for  the  determination  of  eigenvalues.*  We  will  show  here  that 
the  nature  of  the  iterative  procedure  for  evaluating  the  eigenvalues  is 
actually  a  minimization  of  the  Rayleigh  quotient. 

The  equations  of  equilibrium  solved  in  each  step  are  given  by 


(2) 


where  {u,}  represents  all  displacements  other  than  those  on  the  outer 
boundary  Ry,.  During  the  iterative  process,  we  solve 

{u?)  =  -[K,]-'[K,J{uiJ  (3) 

where 

By  further  partitioning  ofeqn.  (2),  {U(}^  =  Lu„,  u„,,  u„J  it  is  possible*  to  write 
eqn.  (3)  as 

[»]{«•"«.}  = (5) 
where  [/I]  and  [5]  are  square  matrices.  Therefore,  eqns  (4)  and  (5)  form 
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a  statement  of  the  generalized  eigenvalue  problem.^  Inverting  the  matrix 
[-4],  we  can  write  eqns  (4)  and  (5)  as 

=  (6) 

where  the  eigenvalue  A  can  be  identified  with  the  radial  function  part  of  eqn. 
(lb),  i.e.  Therefore, 

(7) 

The  displacements  {u^  J  or  eigenfunctions  of  eqn.  (6),  are  identified  with  the 
angular  distribution  of  the  asymptotic  field,  thus 

{u,.}  =  mC,)  (8) 

The  matrix  [T]  is  referred  to  as  the  transfer  matrix.®  It  is  possible  to  show 
that  when  the  matrix  [T]  is  symmetric,  the  singularity  is  real  and  when  it  is 
non-symmetric,  the  singularity  3  is  complex. 

Although  we  have  identified  the  eigenvalue  with  a  power  type  singularity 
in  eqn.  (7),  there  is  no  reason  not  to  assume  that 

A=L(r)  (9) 

where  the  asymptotic  field  of  eqn.  ( 1  b)  can  be  written  in  the  separable  form, 

u,(r,0)=L(r)f/(0)  +  L(r)f?(0)  (10) 

The  function  L(r)  could  be  logarithmic  or  power  singularity  or  a  combina¬ 
tion.  As  long  as  the  field  is  separable,  we  can  find  out  the  L(r)  function  by 
fitting  it  to  the  displacements  along  any  radial  line  in  Fig.  1 . 

2.1.  Complex  Eigenvalues 

Complex  eigenvalue  problems  arise  when  the  transfer  matrix  [T]  is 
non-symmetric.  After  a  large  number  of  iterations  (n),  one  can  write  eqn.  (6) 
as 

N 

3 

where  \j,j=  1, . . . ,  iV  are  conjugate  complex  eigenfunctions  of  the  problem. 
They  represent  a  complete  set  N  for  the  matrix  [T]  (order  Ai).  x  ^  x  j  are  the 
corresponding  complex  conjugate  eigenvectors  of  the  first  dominant 
eigenvalue  Aj ,  and  its  conjugate  X, .  a,-  are  conjugate  constants. 

Using  the  Rayleigh  quotient  approach,  it  is  possible  to  show  that  the  last 
summation  in  eqn.  (10),  will  get  small  in  comparison  with  the  first  two  terms 
as  the  number  of  iterations  gets  larger.  Therefore,  after  a  large  number  of 
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iterations,  the  method  will  converge  to  the  first  dominant  mode,  which  is 
represented  by  the  first  two  terms  of  eqn.  (11). 

At  convergence  we  can  also  write  the  following  equation, 

/Si  {“j. } ‘ } + ^ } = 0  (12) 

where  can  be  determined  by  least  squares  method  for  the  over  determined 
system  of  equations.  The  author  determined  by  satisfying  eqn.  (12)  in  the 
sense  of  the  norms,  i.e.  dot  products  are  used.  These  equations  were  solved 
to  calculate  ,  ^2  ^^d  ^3,  each  equation  was  obtained  by  multiplying  eqn. 
(12)  by  one  of  the  iteration  vectors. 

Using  eqns  (11)  and  (12),  we  get  the  characteristic  equation  for  A,  given 
by: 

/Si+/S2A,  +  ^3^?  =  0  (13) 

from  which  A,  can  be  obtained  and  the  singularity  can  be  calculated  from 

eqn.  (7)  or  (9).  For  a  power  singularity 

=  +  (14) 

or 

i*i  +  i»ji=(Rb/R,)“[cos  (sln(Rb/R,))  +  isin(eln(Rb/R,))]  (15) 

from  which  the  real  and  imaginary  parts,  a  and  e,  can  be  evaluated. 

The  dominant  eigenvector  can  be  calculated  from  two  consecutive 
iterations,  (n)  and  («+  1).  Therefore, 

*1  =  '/i  K,}  +  {uj.}  -/c„+  ,{ui^ ' })  (16) 

where  /c„+i  is  the  scaling  factor  during  that  iteration. 

2.2.  Global  Interpretation  of  FEIM 

The  FEIM  can  also  be  interpreted  as  a  substructuring  method.  Each 
iteration  represents  a  substructure  at  the  radius  R,.  This  means  that  the 
iterations  are  continuously  zeroing-in  on  the  singular  point  (e.g.  crack  tip) 
in  every  successive  analysis.  It  can  be  shown  that  the  stiffness  from  iteration 
(n)  to  iteration  («  + 1)  should  be  scaled  in  the  following  manner, 

[X]"^‘=(R,/Rb)^[A]"  (17) 

Therefore,  if  the  initial  displacements  were  obtained  from  those  of 
a  model  representing  a  structure  or  test  specimen,  the  FEIM,  in  addition  to 
evaluating  the  singularity  field,  will  also  give  the  stress  intensity  factors.  Of 
course,  proper  multiplications  of  the  results  using  eqn.  (17)  will  be  needed. 
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For  the  case  of  complex  eigenvalues,  the  stress  intensities  are  calculated 
from  eqns  (16)  and  (1),  therefore 


x  =  KFJfi)  and  ct  =  x,y  (18a) 

where  the  stress  intensity  K,  in  this  case  is  complex.  It  can  be  written  as; 

K^iK^  +  iK^)  and  F,=f+ig  (18b) 

If  we  also  write  x  in  its  real  and  imaginary  parts,  we  get, 

x=z+iw=(X,+i/:2)(f+ig)  (19) 

Using  dot  products  we  get 

/Ci=(z-f+wg)/(f-f+g-g)  (20a) 

^2  =  -(*‘g  +  w-f)/(f-f+g-g)  (20b) 

For  real  eigenvalues  the  above  equations  reduce  to 


/C,=(z-f)/(f-f) 


(21) 


where  f  is  the  analytical  expression  for  the  functional  variation  in  6  for  the 
mode  of  fracture  in  interest  (I,  II,  or  III).  Equations  (20)  and  (21)  represent 
a  powerful  method  for  evaluating  stress  intensities  for  complex  and  real 
asymptotic  fields.  It  should  be  noted  that  its  accuracy  is  like  integral 
methods. 


2.3.  Verification  of  FEIM 

The  FEIM  was  verified  for  cracks  in  anisotropic  materials.^  Tests  on 
singularities  at  interfaces  included  cracks  perpendicular  to  interfaces,*  free 
surface  meeting  an  interface,*  and  crack  along  an  interface  of  dissimilar 
media.* 

The  results  in  all  these  cases,  where  analytical  solutions  were  available, 
were  extremely  accurate. 

The  most  severe  test*  was  that  of  the  crack  along  the  interface  of 
dissimilar  media.®"  In  that  case,  the  real  and  imaginary  parts  of  eqn.  (14) 
for  the  power  singularity  were  given  by 

a+i£=0-500  004  3-i(0075  78603)  (22) 

The  analytical  solution‘s  for  dissimilar  materials  with  properties 
G,/Gj  =  l/10,  v,  =  v2  =  0-3 
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is  given  by 


a  +  ie=0-5-i(0075  81178) 
where  e  is  calculated  from  the  equation, 


1  r(3-4vi)/Gi  +  l/G/ 

2it  (3— 4v2)/G2  +  1/Gj 


(23) 

(24) 


2.4.  Applications 

The  Finite  Element  Iterative  Method  (FEIM)  as  outlined  here  exhibits 
all  of  the  generalities  of  the  finite  element  method.  Thus  problems  involving 
two-  and  three-dimensional  singularities  as  well  as  various  material 
anisotropy  are  easily  accommodated.  Problems  of  singularity  at  interfaces 
of  isotropic  as  well  as  anisotropic  materials  can  be  handled  easily  with  the 
method  using  the  complex  eigenvalue  method  for  the  determination  of  the 
asymptotic  field.  It  was  shown  that  as  long  as  the  field  can  be  written  in 
a  separable  form  for  the  r  and  6  functions,  the  method  will  work. 

Solutions  of  several  problems  in  composite  materials  are  underway. 
These  include  singularities  at  delaminations,* ^  free  edge  singularities'^ 
at  cut-outs  and  lap  joints,  and  intersection  of  a  crack  between  dissimilar 
media  with  the  free  surface.  The  latter  problem  was  solved  for  homo¬ 
geneous  isotropic  media. '*•'*  However,  some  of  these  results  are  not  in 
agreement  with  the  rest. 


3.  CONCLUSIONS 

It  was  shown  that  the  finite  element  iterative  method  is  a  statement  of  the 
eigenvalue  problem  of  the  asymptotic  singular  field.  The  FEIM  uses 
existing  general  purpose  Finite  Element  Programs  with  little  or  no 
modification,  and  therefore  represents  a  powerful  tool  for  the  practical 
engineer.  The  method  has  all  the  generalities  of  the  finite  element  method 
for  2-D,  3-D,  material  anisotropy,  and  orientation. 

We  have  also  shown  that  the  method  is  also  a  global  method,  and 
whence  the  stress  intensities  for  real  as  well  as  complex  fields  can  be 
calculated. 
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ABSTRACT 

The  transition  from  slow  crack  growth  (ductile  behaviour)  to  rapid  crack 
propagation  (brittle  behaviour)  by  varying  the  structural  size  scale  and 
keeping  the  structural  shape  constant,  is  explained  in  terms  of  dimensional 
analysis  and  described  with  the  application  of  the  Strain  Energy  Density 
Theory.  On  the  other  hand,  in  spite  of  the  apparent  variability  of  the 
constitutive  behaviour  as  a  function  of  size  of  material  element,  the  damage 
mechanics  is  proposed  as  an  invariant  process  at  the  microscale,  which  can  be 
revealed  experimentally  by  temperature  fluctuations. 


1.  INTRODUCTION 

In  the  last  few  years  a  great  effort  has  been  made  by  researchers  and 
institutions  to  explain  two  peculiar  and  recurrent  phenomena  in  material 
strength: 

(1)  size  effect; 

(2)  stable  crack  growth. 

As  regards  the  former,  if  the  structural  size  scale  varies,  the  structural  shape 
being  constant,  the  mechanical  behaviour  of  the  structure  decidedly 
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changes  from  the  very  brittle  to  the  very  ductile  (Fig.  1).  This  effect  can  be 
explained  only  by  the  application  of  physical  similitude  and  scale  modelling 
concepts.'  As  regards  the  latter  phenomenon,  it  represents  only  a  local 
instability  at  the  crack  tip  and  its  nature  is  totally  different  from  that  of 
unstable  crack  propagation.® 


STRUCTURAL 

B£HAVIOUR 

CRACK  GROWTH 
PROCESS 

© 

BRITTLE 

UNSTABLE 

© 

ductile-brittle 

STABLE-UNSTABLE 

© 

DUCTILE 

STABLE 

Fig,  1.  Size-scale  transition  from  brittle  fracture  to  plastic  collapse. 


Stable  crack  growth  may  occur  both  under  monotonic  or  repeated 
loading  and  may  precede  or  follow  the  unstable  crack  propagation.  The 
fundamental  laws  governing  the  transition  from  slow  to  rapid  crack 
propagation,  and  vice  versa,  should  be  very  general  and  applicable  to  very 
simple  as  well  as  to  very  complex  structures,  so  that  it  will  be  easy  and 
consistent  to  extrapolate  the  results  obtained  from  small  specimens  to  the 
project  of  large  structures. 

The  size-scale  transition  from  plastic  collapse  to  brittle  fracture  will  be 
examined  on  the  basis  of  dimensional  analysis  and  described  by  the 
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application  of  the  Strain  Energy  Density  Theory."'  ®.  Then,  the  size  effects 
on  material  strength,  toughness  and  ductility  will  be  discussed  through  the 
scale-invariance  of  the  non-dimensional  crack  growth  resistance  curves. 
Eventually,  the  cooling-heating  effect  in  a  heterogeneous  material  subjected 
to  repeated  loading  will  be  introduced  and  explained  as  a  transition  from 
order  to  disorder  (or  damage)  according  to  the  thermodynamics  of 
irreversible  processes. 


2.  DIMENSIONAL  TRANSITION  FROM  PLASTIC  COLLAPSE 
TO  BRITTLE  FRACTURE 

Two  fundamental  questions  arise  dramatically. 

(1)  Are  the  data  coming  from  small-scale  specimens  related  to  the 
collapse  conditions  in  large-scale  structures? 

(2)  If  the  ductile  fracture  in  small-scale  specimens  is  not  completely 
obscured  by  the  plastic  flow  collapse  at  the  ligament,  how  is  it 
possible  to  put  the  former  in  connection  with  the  brittle  fracture  in 
large-scale  structures? 

The  competition  between  collapses  of  a  different  nature  can  be  empha¬ 
sized  with  the  application  of  dimensional  analysis  and  considering  the 
maximum  loads  derived  from  LEFM  and  limit  analysis  respectively.  The 
transition  from  ductile  to  brittle  behaviour  is  governed  by  a  non-dimen¬ 
sional  brittleness  number  which  is  a  function  of  material  properties  and 
structure  size  scale.  A  true  separation  collapse  occurs  only  with  relatively 
low  fracture  toughness  and/or  large  structure  size. 

Due  to  the  different  physical  dimensions  of  yield  strength,  and 
fracture  toughness,  scale  effects  are  always  present  in  the  usual  fracture 
testing  of  common  engineering  materials.  This  means  that,  for  the  usual  size 
scale  of  the  laboratory  specimen,  the  plastic  collapse  at  the  ligament  tends 
to  anticipate  and  obscure  the  brittle  crack  propagation.  Such  a  competition 
between  different  types  of  collapses  can  easily  be  shown  by  considering  the 
expression  for  the  stress-intensity  factor  in  a  centre  cracked  slab  (Fig.  2): 


where  the  shape  function  /  is; 


(1) 
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Fig.  2.  Interaction  between  brittle  crack  propagation  and  plastic  collapse,  by  varying  the 
brittleness  number  s  =  K,c/(T,y2b. 


At  the  crack  propagation  condition  eqn.  (1)  becomes; 


f<ic~  yj 


If  both  members  of  eqn.  (2)  are  divided  by  o^y/zb,  we  obtain; 


!^=s=^ 


max  fl  ^ 

Oy  V  2h-^\f> 


where  s  is  a  dimensionless  number  able  to  describe  the  brittleness  of  the 
system  and  where  both  material  properties  and  specimen  size  appear. 
Rearranging  eqn.  (3)  gives: 

/  na„\ 


2b 


(4) 
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On  the  other  hand,  it  is  possible  to  consider  even  the  non-dimensional  load 
producing  plastic  collapse  at  the  ligament  (b—aj: 


Equations  (4)  and  (5)  are  plotted  in  Fig.  2  as  functions  of  the  crack  depth 
a^/b.  While  the  former  provides  a  family  of  curves  by  varying  the  brittleness 
number  s,  the  latter  is  represented  by  a  unique  curve  (thick  line).  It  is  easy  to 
realize  that  plastic  collapse  at  the  ligament  precedes  crack  propagation  for 
each  initial  crack  depth  when  the  brittleness  number  s  is  higher  than  the 
critical  value  =  0-54.  For  lower  s  numbers,  plastic  collapse  anticipates 
crack  propagation  only  for  initial  crack  depths  external  to  a  certain 
interval.  This  means  that  a  real  separation  phenomenon  occurs  only  with 
a  relatively  low  fracture  toughness,  high  yield  strength  and/or  large 
structure  size.  Not  the  single  values  of  K^,  and  b,  but  only  their  function 
s  —  see  eqn.  (3)  —  determines  the  nature  of  the  collapse  mechanism. 

3.  STRAIN  ENERGY  DENSITY  THEORY 

3.1.  Mechanical  Damage  and  Strain-Softening  Behaviour 

Damage  of  the  material  at  the  crack  tip  and  crack  growth  increments  will 
be  computed  on  the  basis  of  a  uniaxial  bilinear  elastic-softening  stress-strain 
relation  (Fig.  3(a)).  If  the  loading  is  relaxed  when  the  representative  point 
is  in  A,  the  unloading  is  assumed  to  occur  along  the  line  AO,  so  that  the  new 
bilinear  constitutive  relation  is  the  line  OAF.  No  permanent  deformation  is 
allowed  by  such  a  model,  but  only  the  degradation  of  the  elastic  modulus. 
The  present  model  simulates  the  mechanical  damage  by  decreasing  elastic 
modulus,  £,  and  strain  energy  density  which  can  be  absorbed  by  a  material 
element.  In  fact,  while  for  a  non-damaged  material  element  the  critical  value 
of  the  strain  energy  density,  (dH^dF)^,  is  equal  to  the  area  OUF  (Fig.  3(a)) 
for  a  damaged  material  element  with  representative  point  in  A,  the 
decreased  critical  value,  (d  H^d  V)*,  is  equal  to  the  area  OAF.  In  addition,  as 
is  shown  in  Fig.  3(a),  the  area  OUA  represents  the  dissipated  strain  energy 
density,  (dW/dV)^,  OAB  the  recoverable  strain  energy  density,  (dW/dF)^, 
and  BAF  the  additional  strain  energy  density,  {dW/dV\. 

The  described  model  will  be  extended  to  the  three-dimensional  stress 
conditions,  using  the  current  value  of  the  absorbed  strain  energy  density, 
(dW/dV),  as  a  measure  of  damage.  In  other  words,  the  effective  elastic 
modulus,  £♦,  and  the  decreased  critical  value  of  strain  energy  density, 
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Fig.  3.  Strain-softening  constitutive  law  of  a  concrete-like  material  in  tension:  (a)  assumed 
bilinear  relation;  (b)  numerical  damage  simulation. 


(dWI6V)*,  will  be  considered  as  functions  of  the  absorbed  strain  energy 
density,  (AW/^V),  being; 
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Such  functions  in  the  uniaxial  case  are: 


/dlV\ 

OUAB-area  -♦  AO-slope,  i.e.  I  E*  (6a) 

/dW\  /diV\* 

OUAB-area  -►  OAF-area,  i.e.  ( 

Stress  and  strain  in  the  softening  condition  A  (Fig.  3(a)),  can  be  expressed 
in  terms  of  stress  and  strain  in  the  ultimate  and  fracture  conditions: 


a  =  E*£= - -  (7) 

(£f-£„)  +  |^ 

Through  eqns.  (6)  and  (7)  it  is  simple  to  express  the  absorbed  strain  energy 
density,  (d  Wjd  V),  and  the  decreased  critical  value,  (d  Wjd  K)* ,  as  functions  of 
the  constitutive  parameters,  o-y,  £„,  £f ,  and  of  the  effective  Young’s  modulus, 
£*: 


(dW\  1  ^ 

(dW\  1 

[dvl-2''' 


fdW\  /'dW\  /'dW\  1 

fdW\*  fdW\  fdW\  1 


e  +  <rej 


The  relations  (8)  will  be  discretized  using  25  different  values  of  the  elastic 
modulus: 

£*(n)=^^^£,  for  n  =  l,2,...,25  (9) 


In  Fig.  3(b),  the  results  of  different  tensile  test  numerical  simulations  are 
reported.  Two  strain-controlled  loading  processes  are  carried  out  and  the 
strain  increment  in  one  case  is  twice  as  in  the  other  one.  The  discretized 
stress-strain  relations  are  displayed  and  the  comparison  with  the  assumed 
a-e  constitutive  law  is  also  shown.  It  can  be  proved  that,  when  the  effective 
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elastic  modulus,  £*,  varies  continuously  and  the  strain  increment  Ae  tends 
to  zero,  the  assumed  bilinear  (t—e  variation  (dashed  line  of  Fig.  3(b))  is 
exactly  reproduced  by  the  numerical  damage  simulation. 


3.2.  Slow  Crack  Growth  versus  Unstable  Crack  Propagation 

In  order  to  evaluate  the  crack  growth  increment  at  each  loading  step,  the 
Strain  Energy  Density  Theory  will  be  applied,  as  proposed  by  Sih.’  *  It  is 
based  on  the  following  fundamental  assumptions. 

(1)  The  stress  field  in  the  vicinity  of  the  crack  tip  cannot  be  described  in 
analytical  terms  because  of  the  relative  heterogeneity  of  the  material. 
A  minimum  distance,  r„,  does  exist  below  which  it  is  a  nonsense  to 
study  the  mechanical  behaviour  of  the  material  from  a  ‘continuum 
mechanics’  point  of  view  and  to  consider  macroscopic  crack  growth 
increments. 

(2)  Out  of  such  a  core  region  of  radius  r^,  the  strain  energy  density  field 
can  always  be  described  by  means  of  the  following  general  relation¬ 
ship: 


where  the  strain  energy  density  factor,  S,  is  generally  a  function  of  the 
three  space  coordinates. 

(3)  According  to  Beltrami’s  criterion,  all  the  material  elements  in  front  of 
the  crack  tip,  where  the  strain  energy  density  is  higher  than  the 
critical  value,  (dW^dK)*,  fail. 

(4)  When  the  following  condition  holds: 

^a  =  r*  =  SJm|<^V)^  (11) 

the  crack  may  be  considered  as  arrested,  at  least  from  a  macroscop- 
ical  point  of  view.  On  the  other  hand,  when  the  crack  growth 
increment  is: 

Aa  =  r?  =  S,/(dlT/dK)?  (12) 

the  unstable  crack  propagation  takes  place.  is  a  material  constant 
and  represents  the  strength  of  the  material  against  rapid  and 
uncontrollable  crack  propagation.  is  connected  with  the  critical 
value  of  the  stress-intensity  factor,  through  the  following 
equation  (plane  strain  condition).^ 


S  = 


(I-hv)(l-2v) 


K 


2 

1C 


2kE 


(13) 
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33.  Centre  Cracked  Slab  in  Tension 

A  centre  cracked  slab  in  tension  (Fig.  4(a))  is  analysed  by  using  the 
Axisymmetric/Planar  Elastic  Structures  (APES)  finite  element  program.’  It 
is  a  computer  program  which  incorporates  12-noded  quadrilateral  isopara¬ 
metric  elements  allowing  for  cubic  displacement  fields  and  quadratic  stress 
and  strain  fields  within  each  element.  The  strain  energy  density 
singularity  in  the  vicinity  of  the  crack  tip  is  embedded  in  the  solution 
through  the  use  of  nodal  spacing  on  the  element  sides  adjacent  to  the 
crack  tip.  The  idealization  of  Fig.  4(b)  utilizes  309  nodes  and  52  elements 
and  is  considered  in  a  condition  of  plane  strain. 


niD” 


T 


II 

g 


Fig.  4.  Centre  cracked  slab  in  tension. 


The  strain  energy  density  criterion  is  applied  to  the  tension  test  specimen, 
considering  a  strain-controlled  loading  process.  The  stress-strain  responses 
for  three  different  initial  crack  lengths  are  displayed  in  Fig.  5.  The  load 
carrying  capacity  decreases  by  increasing  the  initial  crack  length. 

In  Fig.  6,  the  stress  a  is  represented  against  the  crack  growth,  2(a  —  a,).  At 
the  first  steps,  the  stress  increases  while  the  crack  grows.  Then,  after 
reaching  a  maximum,  the  stress  decreases  and  attains  the  value  zero  when 
the  whole  ligament  is  separated.  The  maximum  represents  the  transition 
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Fig.  6.  Stress  versus  crack  growth  plots  for  three  different  initial  crack  lengths. 


between  stable  and  unstable  structural  behaviour.  On  the  other  hand,  the 
transition  between  stable  and  unstable  crack  propagation  depends  on  the 
achievement  of  the  critical  value  of  the  strain  energy  density  factor,  .  Such 
a  value  has  physical  dimensions  different  from  those  of  (AWIdV)^,  and 
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produces  the  scale  effects  recurrent  in  Fracture  Mechanics.  In  the  reality, 
the  crack  instability  may  precede  or  follow  the  structural  instability.  This 
mainly  depends  on  the  structural  size  scale. 

Because  material  damage  and  crack  growth  occur  in  a  non-self-similar 
fashion  for  each  step  of  loading,  specimens  of  different  sizes  appear  to 
behave  differently.  The  application  of  Buckingham's  theorem  for  physical 
similitude  and  scale  modelling  gives; 

rr  ^  '  flo 

_  [dvl 

where  material  toughness,  (dWjAV)^,  and  specimen  width,  b,  have  been  used 
as  the  fundamental  quantities.  The  stress,  a,  may  be  regarded  as  a  function 
of  the  strain  e  only,  if  all  other  ratios  are  kept  constant. 

In  the  same  way,  it  is  possible  to  define  a  dimensionless  strain  energy 
density  factor; 


V !  I  ^ 

'  b'b'  b 


(15) 


Function  S  can  be  regarded  as  linear  in  a/b^ 

S  dS/da  a-a„  S„ 

TdwYT  ~  (dIF/dnc  ^  TdTr V 

which  may  obviously  be  rearranged  into  the  form; 


The  constants  A  and  B  are  dimensionless  and  scale  independent.  It  follows 
that  the  slope  of  the  S-a  diagram  is  constant  varying  the  scale  and  the 
intercept,  S„,  is  proportional  to  the  scale  b. 

Figure  7  shows  the  straight  line  plots  of  S  versus  crack  growth 
for  increasing  size  b{a„/b  =  0-3).  The  critical  crack  growth  decreases 
with  increasing  specimen  size.  For  example,  with  the  critical  value  = 
8x10"^  kg/cm,  the  limiting  size  is  6  =  240  cm.  Beyond  this  size,  stable  crack 
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growth  ceases  to  occur  and  failure  corresponds  to  unstable  crack  propaga¬ 
tion  or  catastrophic  fracture. 

Figure  8  presents  the  relations  between  stress  and  strain  (a,/h=0-3).The 
vertical  lines  with  arrows  indicate  the  limiting  values  of  e  as  the  critical 


CRACK  GROWTH.  2s  -  2Sg  (cm) 

Fro.  7,  Strain  energy  density  factor  versus  crack  growth  plots  by  varying  size  b  (a  Jb =()•%). 
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Fig,  8.  Stress  versus  strain  relations  by  varying  size  b  (a„/b=0-3). 
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Strain  energy  density  factor,  =  8  x  10“  ^  kg/cm,  is  reached.  Crack  instabil¬ 
ity  occurs  for  smaller  strains  as  the  size  b  increases.  This  is  obvious,  since  the 
initial  crack  length  also  increases,  the  ratio  a^/b =0-3  being  constant.  The 
structural  instability  occurs  before  crack  instability  only  for  b  ^  80  cm  (Fig. 
8).  When,  for  example,  b  =  120  cm,  softening  behaviour  is  not  present  and 
the  crack  starts  to  spread  in  an  unstable  manner,  stress  a  still  being  in  the 
ascending  stage. 

It  is  also  interesting  to  consider  the  maximum  stress  resulting  from 
the  Linear  Elastic  Fracture  Mechanics  (LEFM)  solution,  eqn.  (4),  and 
coming  from  the  limit  analysis  at  the  ligament,  eqn.  (5).  Normalizing 
the  strengths  ff^J^and  obtained  respectively  by  the  present  approach 
and  by  the  limit  analysis,  with  obtained  through  LEFM,  we  can 
evaluate  the  interaction  between  the  different  failure  modes  (Fig.  9).  The 


Fig.  9.  Transition  between  plastic  collapse  and  brittle  crack  propagation  by  varying  size 

b(aJb=0-3}. 

horizontal  line  and  equal  to  100%  represents  the  case 

when  failure  coincides  totally  with  brittle  fracture,  while  the  dashed  line 
represents  the  case  when  failure  coincides  totally  with  plastic  collapse. 
When  the  specimen  size  is  small,  the  simple  formula  in  eqn.  (5)  gives  good 
prediction  based  on  the  ultimate  strength  alone.  On  the  other  hand,  when 
the  specimen  size  is  large,  eqn.  (4)  gives  good  prediction  based  on  the 
stress-intensity  factor  alone.  The  two  extreme  situations  are  then  connected 
by  a  transition.  For  intermediate  sizes,  the  ratio  of  Fig.  9  appears  higher 
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than  one.  This  means  that  a  fictitious  critical  stress-intensity  factor 
K\^  larger  than  the  true  may  be  assumed.  In  this  way,  the  energy¬ 
absorbing  damage  process  at  the  crack  tip  is  taken  into  account. 


4.  SIZE  EFFECTS  ON  STRENGTH, 

TOUGHNESS  AND  DUCTILITY 

Equation  (16)  may  obviously  be  rearranged  into  the  form: 

S=S'{^-U  +  So  (18) 

the  constants  S'  and  being  dimensionless  and  scale  independent.  On  the 
other  hand,  it  is  apparent  that  the  quantity  S^I('AWIAV)^b  must  also  enter 
into  the  dimensional  analysis  in  eqn.  (14).  In  fact,  for  estimating  it 
suffices  to  consider; 


=  n{S*) 


(19) 


in  which  S*  is  a  brittleness  number  analogous  to  that  defined  in  eqn.  (3): 

S. 


S*  = 


i;c 

b 


(20) 


Hence,  all  geometrically  similar  structures  can  be  regarded  as  governed  by 
S*.  This  dimensionless  quantity  can  be  used  to  predict  the  stress  versus 
strain  behaviour  for  all  specimen  sizes.* 

In  conclusion  and  recalling  eqn.  (18),  it  is  possible  to  state  that  above  the 
size  (Fig.  10): 


^max 


s„ 


(21) 


stable  crack  growth  ceases  to  occur  and  the  brittle  failure  is  achieved  when 
the  a-E  curve  is  still  in  its  linear  initial  course  (Fig.  8),  whereas  below  the  size 
(Fig.  10): 


(^min 


Sc 


^„-^5'(i-<^j 


dW 


(22) 
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Fig.  10.  Scale  invariance  of  the  non-dimensional  crack  growth  resistance  curve. 


unstable  crack  growth  ceases  to  occur  (even  in  the  softening  stage)  and  the 
progressive  slow  crack  growth  develops  up  to  the  complete  specimen 
separation.  On  the  other  hand,  for  stable  (or  slow)  crack 

growth  is  followed  by  unstable  (or  fast)  crack  propagation.  Ductility 
appears  therefore  as  a  mechanical  property  depending  on  the  size  scale  of 
the  structure,  as  well  as  on  the  fracture  toughness  of  the  material. 


5.  COOLING-HEATING  EFFECT  IN  A  HETEROGENEOUS 
MATERIAL  SUBJECTED  TO  REPEATED  LOADING 

5.1.  Description  of  Material  and  Testing  Procedure 

In  spite  of  the  apparent  variation  of  the  constitutive  behaviour  as 
a  function  of  the  size  of  the  material  element,  the  damage  mechanics  is 
proposed  as  an  invariant  process  at  the  microscale,  which  can  be  revealed 
experimentally  by  temperature  fluctuations. 

The  first  of  a  series  of  tests  on  composite  materials,  such  as  concrete,  is 
described  within  the  framework  of  an  investigation  aimed  at  identifying 
a  thermodynamic  model  for  predicting  the  fatigue  life  of  these  materials 
when  subjected  to  repeated  compression  loading.  The  tests  reveal  a  marked 
‘cooling-heating’  effect  for  very  modest  maximum  compressive  stresses 
=  where is  the  ultimate  compressive  strength).  The  repeated 
loading  cycles  produce  a  decrease  in  the  temperature  of  the  material  at  first. 
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and  then,  as  the  applied  load  increases,  the  temperature  is  seen  to  increase. 

The  behaviour  is  consistent  with  the  findings  of  recent  investigations  on 
steel  and  aluminium  specimens  tested  in  tension.”*"*^  However,  since  we 
are  dealing  with  compression  tests  and  repeated  loading,  a  further 
investigation  was  considered  useful. 

It  is  verified  that  the  ‘cooling-heating’  effect  depends  on  the  value  of  the 
maximum  applied  stress  The  test  specimens  subjected  to  stress 

present  the  heating  stage  only  (in  this  case  the  cooling  phase 
may  occur  and  yet  be  concealed).  In  the  experimental  program,  carried  out 
at  the  Materials  Testing  Laboratory  of  the  Structural  Engineering 
Department  of  the  Politecnico  di  Torino,  another  material  was  also 
considered  (plexiglass)  for  comparison  purposes. 

The  first  material  is  concrete  made  with  350  kg/m^  Portland  cement  and 
alluvial  aggregate  with  max.  diameter  (^  =  8  mm  =  0-3 1  in.  The  compressive 
strength  was  found  to  bef^  =  37  MPa.  The  test  is  performed  on  a  cylindrical 
specimen  (H  =  D  =  2  i  cm  =  l  lOin),  thermally  insulated  from  the  environ¬ 
ment  by  means  of  polystyrene  applied  to  end  and  side  faces  (Fig.  1 1).  The 
temperature  of  the  specimen  is  measured  by  means  of  two  thermocouples, 
Tj  and  T2,  applied  to  the  side  faces,  and  a  third  thermocouple  T3  is  used  to 
determine  the  temperature  of  the  environment. 


Fig.  11.  Concrete  specimen  under  compression  loading  cycles.  Testing  apparatus. 
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The  comparative  specimen  of  plexiglass  consists  of  a  prism  (3x3x5  cm) 
tested  under  identical  conditions. 

The  temperature  versus  time  curves  are  plotted,  the  specimens  being 
subjected  to  rejieated  compression  loading  by  means  of  a  MTS  machine 
(max.  loading  capacity  of  1(X)  000  N)  and  by  measuring  the  temperature 
values  by  means  of  a  Compact  Logger  3440  with  a  sensitivity  of  01  K. 

5.2.  Results  and  Comments 

Table  1  reports  the  concrete  specimen  temperature  values  in  relation  to 
load  variations.  The  load  is  characterized  by  the  and  values  and 
the  oscillation  frequency. 

The  temperature  variation  Ar(AT=7J  — 7^) = specimen  temperature  — 
environmental  temperature=((T,  +  T2)/2)  — Tj  is  plotted  in  Fig.  12. 

For  load  I(ffm„//c  =  0- 1 35)  at  an  oscillation  frequency  of  2  Hz  we  observe 
a  marked  drop  in  temperature  (cooling  stage)  which  becomes  stable 
around  an  asymptotic  value  AT s;  1  K.  This  value  remains  unchanged  with 
load  II  (ffm„//c =0-27).  A  reversal  in  temperature  is  seen  to  occur  when  load 


Fig.  12.  Temperature  versus  time  diagram  for  a  concrete  specimen  under  compression 
loading  cycles  at  frequencies  =  2-6  Hz. 


Transfer  of  Small  Specimen  Data  to  Full-Size  Components 
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III  is  applied  At  this  last  stage,  the  effects  of  the  load 

oscillation  frequency  on  temperature  values  can  be  observed.  The  tem¬ 
perature,  after  reaching  a  plateau,  increases  when  the  frequency  is  raised 
from  2  to  6  Hz. 

Figure  13  shows  the  results  obtained  for  the  plexiglass  specimen.  The 


Fig.  13.  Temperature  versus  time  diagram  fora  PMMA  specimen  under  compression  loading 
cycles  at  a  frequency  of  2  Hz. 

same  testing  procedure  adopted  for  the  concrete  test  was  used.  The  cooling 
stage  does  not  appear,  at  least  according  to  the  sensitivity  of  the  equipment 
used  (±01  K).  Beginning  with  load  III  ((7„,„ = 4-4  M  Pa)  a  marked  increase 
in  temperature  occurs. 

It  should  be  observed  that  the  stress-strain  curves  obtained  for  plexiglass 
during  the  test  by  means  of  electrical  strain  gauges  turn  out  to  be  perfectly 
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overlapping  with  increasing  cycle  numbers.  This  means  that  energy 
dissipation  can  be  neglected.  The  opposite  happens  in  the  case  of  c  .^ncrete, 
for  which  the  hysteretic  curve  shifts  in  the  direction  of  positive  strains,  with 
permanent  deformation. 

The  ‘cooling-heating’  effect  may  be  explained  in  terms  of  the  concepts  of 
order  and  disorder  formulated  in  the  ‘thermodynamics  of  irreversible 
processes’.  The  material’s  cooling  coincides  with  the  period  in  which  order 
is  established  in  the  system  with  respect  to  the  initial  state.  Physically,  this 
could  correspond  to  an  increase  in  the  material’s  homogeneity.  In  this  way, 
the  differences  between  concrete  and  plexiglass  may  be  evidenced.  In 
plexiglass,  in  fact,  the  homogenisation  stage  is  much  less  pronounced  than 
in  concrete:  therefore,  cooling  is  not  detected.  The  temperature  versus  time 
curve  represents  a  possibility  of  measuring  the  state  of  order  and  disorder 
and  therefore  it  becomes  a  precious  tool  to  forecast  the  fatigue  life  of 
stiuctural  materials.  The  minimum  in  the  temperature  versus  time  diagram 
seems  to  represent  the  transition  between  order  and  disorder,  from 
a  thermodynamic  point  of  view,  or  between  reversible  deformation  and 
irreversible  damage,  from  a  mechanical  point  of  view. 
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ABSTRACT 

Herein  is  contained  a  theoretical  assessment  of  the  resulting  damage  in 
composite  laminates  subjected  to  low-velocity  impact.  The  composite  is 
considered  to  be  planar  isotropic  as  found  in  Chopped  Strand  Mat  (  CSM  ) 
and  Sheet  Moulding  Compounds  ( SMC ).  The  proposed  theoretical  method¬ 
ology  will,  although  approximate  by  necessity,  be  shown  to  reasonably  predict 
the  loss  in  strength  and  stiffness  which  occurs  when  such  composites  are 
subjected  to  impact  damage.  Reasonable  comparison  with  experimental 
studies  on  glass/ polyester  CSM  samples  will  be  demonstrated. 


NOTATION 


W, 

% 

Vv 

K, 

F(t) 

M 

Fi,F2,etc. 


T 


Displacement  of  spherical  impactor 
Displacement  of  test  plate 
Impacting  velocity 
Volume  of  damage 
Force/time  relationship 
Indentor  mass 
Material  failure  stresses 
Peak  force  during  impact 
Available  energy  to  cause  damage 
Plate  kinetic  energy 
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u 

t 

X 

m,  n,  i,j 
a,  b 
c 

<7 

w 

h 

X,  y,z 
r,  6,  z 
t/r,  t/, 

P 

a> 

a 

*0 

a*,b* 

<^US 

<^0 


Plate  flexural  energy 
Maximum  time  being  considered 
Any  typical  intermediate  time,  i.e.  0  ^  r  ^  t 
Integer  counters 

Test  specimen  planform  dimensions 
Planform  radius  of  indentation 
Power  index 

Lateral  plate  displacement 

Plate  thickness 

Plate  cartesian  co-ordinates 

Polar  co-ordinates  in  damaged  region 

Local  displacements  at  impacted  region  in  the  radial  and 

transverse  directions  respectively 

Mass  density  of  composite 

Vibration  frequency 

Typical  plate  indentation 

Maximum  plate  indentation 

Permanent  plate  indentation 

Characteristic  crack  dimensions 

Damaged  specimen  residual  strength 

Undamaged  specimen  residual  strength 

Peak  contact  stress 

Transverse  or  through-thickness  stress 


1.  INTRODUCTION 

At  the  outset  it  should  be  stated  that  presently  there  exists  no  generalised 
means  of  analytically  modelling  the  damage  in  composite  structures 
subjected  to  impact  damage.  This  is  true  for  both  ballistic  and  low-velocity 
impact  damage.  Moreover,  the  generic  term  given  to  studies  of  this  nature, 
damage  tolerance,  is  not  uniquely  defined  and  has  a  number  of  con¬ 
notations  depending  on  the  particular  field  of  study.  For  example,  the  term 
‘damage  tolerance’  has  been  used  to  describe  the  advanced  stages  of 
damage  near  the  fatigue  limit  of  structures  subjected  to  cyclic  loading.  It  has 
also  been  used  to  describe  the  ability  of  structures  to  successfully  carry 
loading  after  localised  impact  damage,  i.e.  their  tolerance  to  damage.  An 
earlier  paper  by  the  authors*  cited  references  reflecting  a  spectrum  of  usage 
of  the  term  damage  tolerance. 

At  this  stage,  it  should  be  appreciated  that  the  effects  of  low-velocity 
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impact  of  composite  materials  will  be  appreciably  different  from  the  more 
traditional  isotropic  equivalents.  An  isotropic  body  will  generally  deform 
elastically,  then  plastically  leaving  permanent  indentation,  whereas  the 
composite  will  undergo  a  multimode  failure,  i.e.  the  impacted  energy  will  be 
transformed  into  matrix  yielding  and  cracking,  delamination,  fibre 
cracking  and  pull-out  or  a  complex  mixture  of  these  failure  modes.  Since 
each  of  these  forms  of  failure  is  difficult  to  analytically  model,  their 
interaction  produces  extreme  mathematical  difficulties.  Thus  the  whole 
concept  of  mathematically  modelling  failure  mechanisms  in  composite 
structures  is  fraught  with  considerable  difficulties  irrespective  of  the 
available  computational  capability.  Hence  the  reason  for  the  present 
workshop  which  critically  addresses  this  problem  and  assesses  the  present 
state  of  the  art  in  this  field. 


2.  THEORY 

The  methodology  proposed  herein  will,  of  necessity,  contain  a  number  of 
simplifying  assumptions.  Although  the  same  general  philosophy  can  be 
applied  to  both  thin  and  thick  structures,  the  former  will  be  specifically 
addressed  as  most  composite  structures  are,  relatively  speaking,  thin. 
Essentially  the  main  difference  in  the  response  of  thin  and  thick  structures 
to  low-velocity  impact  is  flexural  deformations  which  are  prevalent  in  the 
former.  Consequently,  the  energy  consumed  by  such  deformations  must  be 
accounted  for  when  evaluating  the  available  energy  to  cause  damage  of  the 
composite. 

Applying  the  Duhamel  integral  to  the  motion  of  the  impactor,  the  plate 
equations  (1)  and  (2)  are  obtained. 

*^,=  ^0,  -  F(TKt-T)dT  (1) 

^  II  f sin  -  t)  dr  (2) 

m  R  J  0 

Hence  the  indentation  of  the  test  plate  will  be  given  by  (3) 
ac=Vot-^\ F(TKt-T)dT-^yj.—  f  f(T)sincu„,(t-T)dT  (3) 

The  impact  of  a  thick  plate  will  yield  the  force/time  and  displace¬ 
ment/time  relationships  shown  in  Fig.  1,  with  the  equivalent  curves  for 
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Fig.  1.  Typical  force  and  displacement  time  history  for  the  impact  of  a  thick  plate. 


a  thin  plate  given  in  Fig.  2.  Worthy  of  note  is  that  eqn.  (3)  is  only  valid  until 
the  point  of  unloading,  at  which  time  the  load-indentation  law  changes  and 
eqn.  (4),  which  is  the  general  unloading  law,  can  be  used  to  represent  the  left 
hand  side  of  eqn.  (3).  Typical  repeated  loading  and  unloading  of  a  thin 
composite  plate  is  shown  in  Fig.  3. 

rF(T)T'’ 

“  =  (am-«o)  +  ao  (4) 

This  law  was  experimentally  verified  as  diagrammatically  illustrated  in 
Fig.  4.  By  numerically  solving  eqns  (3)  and  (4)  by  the  mean  force  method,  the 
complete  loading  history  can  be  obtained. 

The  energy  available  to  cause  damage  will  not  be  totally  absorbed  by  the 
plate  and  can  be  split  into  three  main  headings;  that  which  causes  damage, 
that  which  induces  plate  vibrations  and  the  rebound  kinetic  energy.  The 
energy  to  cause  damage  is  given  by  eqn.  (5).  The  vibrational  energy  comes 
from  two  parts,  kinetic  and  bending,  eqns  (6)  and  (7),  with  the  total  given 
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Indentation 

Fig.  3.  Loading  and  unloading  of  a  thin  composite  plate. 


K]  =  K]-‘[a,] 
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The  internal  displacements,  eqn.  (11)  can  be  obtained  by  curve  fitting 
a  combination  of  experimental  and  finite  element  solution  points. 


Thus  the  internal  stresses  can  be  evaluated  by  substitution  of  (1 1)  into  (10) 
and  (9). 

The  foregoing  results  are  appropriate  in  the  vicinity  of  the  impacted 
region.  Thereafter,  the  areas  remote  from  damage  can  be  considered  using 
a  cartesian  set  of  co-ordinates. 

The  in-plane  strains  remote  from  the  damaged  region  can  be  written  as: 


d^w 

Y  =-2z  — 
dxdy 

The  plate  stresses  and  strains  can  be  related  by  eqn.  (13). 

=  « 1 1  O’;.  +  +  « 1 

^y  —  ^\2^x^^22^y'^^26^xy 
Vxy  ^16^ X  ^26^y  ^66^ xy 

Therefore,  the  in-plate  stresses  are  given  by  eqn.  (14). 


r  d^w 


d^w 


+  A2 


d^w 


+  2A 


'  ~“^^dxdy\ 


d^w  d^w  d^w 


(12) 


(13) 


(14) 


Theoretical  Modelling  of  Damage  in  Composite  Laminates  141 

From  simple  equilibrium  of  a  typical  plate  element  eqns  (15)  can  be 
formed. 


0Tj,„ 

H - — 

dx 

dy 

dz 

Jjyi. 

dx 

dy 

dz 

Thus  the  corresponding  through-thickness  shear  stresses  are  given  by  eqn. 
(16). 

if  2  1‘^ir  .<  1  i  ^  -1 ..  V  A 


if  2 


d^w  d^w 

r-^  +  (Ai2  +  ^^66)T-2j:  +  '^26 


dx^dy 


Considering  the  specific  case  of  a  rectangular  planform  plate  being 
impacted  by  a  spherical  indentor,  the  total  combined  stresses  can  be  given 
as  eqn.  (17),  with  the  individual  stresses  transformed  into  the  polar 
co-ordinate  system  by  eqns  (18). 


O'TnIiil  —  0^1 


<T,  =  (T,  cos^  0  -I-  (Tj,  sin^  0  -I-  2x^y  sin  9  cos  9 
dg  =  sin^  0  +  <t  cos^  0 -  2t„  sin  0  cos  0 


^ez=-'^xz  sin0-f-T,„  cos  9 
r„=-r^,cos9  +  ty^sm9 
Trj = {Oy  -  ffj  sin  0  cos  0  -1-  T,(j,(cos^  0  -  sin^  0) 

With  the  internal  stresses  evaluated,  the  general  stress  distribution  is  given 
by  Fig.  5,  if  a  suitable  failure  criterion  is  applied  on  three  mutually 
perpendicular  faces  of  an  element  using  eqns  (19). 


o\ 

<Ti<T, 

‘  ^  1 

«^2  ffu 

F? 

F.Fa 

FrF.^, 

(jI 

<^3  023 

F\ 

Fl  F2\ 

O3  ffl3 

Fl 

f.F3^ 

Fl  F^ 
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Using  this  methodology,  a  typical  symmetrical  section  of  the  resulting 
volume  of  damage  will  be  shown  in  Fig.  6.  By  representing  the  characteristic 
dimensions  of  the  volume  of  damage  using  a  semi-elliptical  crack,  as  shown 
in  Fig.  7,  and  making  the  resulting  behaviour  of  the  damaged  and  cracked 
specimens  identical,  a  measure  of  the  residual  properties  of  the  damaged 
specimen  can  be  obtained. 

The  application  of  linear  elastic  fracture  mechanics  gives  the  strain 
energy  release  rate  as  eqn.  (20). 


G,= 


(20) 
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where  the  stress  intensity  factor  is  given  by  eqn.  (21). 

1 


fin  nb*^\ 


nfl*\ 

nb*^\\na*^^^  2h  ) 
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(21) 


Thus  the  residual  strength  can  be  evaluated,  eqn.  (22). 

<^i>s  =  «^us[^I+^^Ke-W]J  (22) 

Providing  a  relationship  between  stress  and  strain  is  available,  the  residual 
stiffness  can  also  be  evaluated. 


3.  RESIDUAL  STRENGTH  INDICATIONS 

Typical  residual  strengths  for  thick  and  thin  plates  are  shown  in  Figs.  8  and 
9,  respectively.  The  latter  figure  shows  a  number  of  experimental  points 
obtained  using  CSM-GRP  specimens  which  corroborate  the  present 


1-0 


Damaged  plate 


10 


Energy  (J) 


Fig.  8.  Residual  damaged  strength  (thick  plates). 
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Fig.  9.  Experimental  and  theoretical  residual  strength  (thin  plates). 


theoretical  methodology.  Also,  the  former  figure  clearly  shown  that 
appreciably  greater  levels  of  damage  energy  are  necessary  to  produce 
significant  loss  of  strength  in  relatively  thick  specimens. 

Clearly  the  present  methodology  has  wide  connotations  with  regard  to 
various  forms  of  composite  materials,  and  indeed,  traditional  materials  of 
construction,  the  only  proviso  being  the  assumption  of  planar  isotropy. 
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ABSTRACT 

A  forthcoming  paper  describes  an  exact  method  of  solution  of  the  three- 
dimensional  elasticity  equations  for  both  the  stretching  and  bending  of 
isotropic  laminated  elastic  plates.  The  theory  is  a  generalisation  to  laminates 
of  the  'exact  plane  stress'  theory for  homogeneous  plates  which  was formulated 
by  Michell  in  1900  and  was  described  by  Love.*  The  essential  feature  is 
that  any  solution  of  the  two-dimensional  classical  elastic  plate  equations  can 
be  used  to  generate  an  exact  three-dimensional  solution.  The  two-dimensional 
solution  may  be  obtained  by  any  of  the  available  methods,  including  numerical 
methods.  Among  other  quantities,  the  interfacial  shear  stress  components  are 
evaluated  exactly. 

Progress  in  extending  this  theory  to  the  more  important  case  of  anisotropic 
laminates  is  described.  In  this  case  exact  closed  form  solutions  cannot  be 
obtained,  but  solutions  in  infinite  series  of powers  of  the  aspect  ratio  have  been 
formulated. 

This  approach  yields  solutions  which  exactly  satisfy  the  three-dimensional 
elasticity  equations  and  lateral  surface  boundary  conditions.  Edge  boundary 
conditions  are  satisfied  only  in  an  average  sense.  However  the  deviation  of  the 
edge  boundary  values  given  by  the  theory  from  those  prescribed  is  a  measure  of 
the  magnitude  of  the  edge  effects  which  may  cause  delamination  and  failure. 

1.  INTRODUCTION 

The  three-dimensional  linear  elastic  stress  distribution  in  a  laminated  plate 
is  basic  to  any  analysis  of  damage  or  failure  in  the  plate.  In  this  paper  we 
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outline  a  procedure  for  determining  such  stress  distributions  for  plates  of 
any  shape  subject  to  any  of  the  standard  edge  boundary  conditions.  Both 
stretching  and  bending  solutions  may  be  dealt  with,  but  for  brevity  we 
concentrate  on  the  problem  of  in-plane  stretching  under  edge  loading. 

With  the  exception  of  a  few  very  special  solutions,  the  existing  theories  for 
elastic  analysis  of  laminated  plates  are  approximate  theories,  and  none  of 
them  appears  to  the  author  to  be  entirely  satisfactory.  A  brief  review  is  given 
in  Ref.  1.  The  most  commonly  used  theory  is  probably  the  one  known  as 
Classical  Laminate  Theory,  which  is  described,  by  e.g.  Christensen.^ 
In  effect,  in  this  theory  the  laminated  plate  is  replaced  by  a  homogeneous 
equivalent  plate  of  the  same  overall  geometry  but  with  elastic  constants 
which  are  suitably  weighted  averages  of  the  elastic  constants  of  the  laminae. 
By  its  nature,  such  a  theory  can  only  yield  averaged  displacement 
components  and  approximate  stress  components.  Whilst  this  knowledge  is 
adequate  for  many  purposes,  for  the  analysis  of  failure  it  is  important  to 
know  the  detailed  through-thickness  distribution  of  stress  and  displace¬ 
ment  so  as  to  determine,  for  example,  the  inter-ply  shear  tractions  which 
have  an  important  bearing  on  the  onset  of  delamination. 

In  a  recent  paper,  Kaprielian  et  a/.‘  have  given  a  three-dimensional 
analysis  for  a  laminate  composed  of  isotropic  laminae.  This  theory  is 
outlined  in  Section  2.  The  outstanding  feature  of  the  theory  is  that  the 
three-dimensional  solution  for  the  laminate  is  generated,  in  a  very  simple 
manner,  by  the  two-dimensional  classical  laminate  theory  solution  for  the 
equivalent  plate. 

The  case  of  greater  practical  interest  is  that  in  which  the  laminae  are 
anisotropic  and  differently  oriented.  The  theory  of  Section  2  does  not  extend 
directly  to  the  case  of  anisotropic  laminae,  but  it  can  be  adapted  to  this  case, 
to  any  required  degree  of  accuracy.  Work  on  the  problem  of  anisotropic 
laminae  is  still  incomplete,  but  a  preliminary  report  on  it  is  given  in  Section 
3.  It  is  intended  that  a  full  account  will  be  published  elsewhere. 

The  only  restriction  on  these  theories  is  that,  as  in  all  plate  theories,  edge 
boundary  conditions  can  be  satisfied  only  in  an  average  fashion,  rather  than 
point  by  point.  A  brief  discussion  of  edge  effects  is  presented  in  Section  4. 


2.  ISOTROPIC  LAMINATES 

This  section  gives  a  brief  summary  of  the  theory  which  is  described  in  detail 
in  Ref  1.  We  first  state  some  results  of  classical  isotropic  thin  plate  theory. 
Our  plate  is  supposed  to  lie  in  the  X,  Y  plane  of  a  system  of  rectangular 
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Cartesian  coordinates  X,  Y,  Z.  Relative  to  this  system,  displacement 

components  are  denoted  by  U,  V,  IK  and  stress  components  by 

a„.  Then  for  in-plate  stretching  deformations  the  governing  equations  are 

2(»/+l)A,;f-n,y=0,  2(^-l-l)A^-Hn.;^=0  (2.1) 

where 

A  =  C7_^+r^,  ft=F;f-£7y  (2.2) 

and 

i/  =  A/(2-f-2A/)=v/(l-v)  (2.3) 

where  C7,  V  are  average  in-plane  displacements,  A,  fi  are  the  Lame  elastic 
constants,  v  is  Poisson’s  ratio,  and  subscripts  following  commas  denote 
partial  derivatives.  It  follows  from  (2.1)  that 

V^A=0,  V2Q=0  (2.4) 

where  is  the  two-dimensional  Laplacian  operator.  For  bending 
deformations  under  edge  loading,  the  classical  result  is 

(2.5) 

where  W  is  the  deflection  of  the  mid-plane.  When  V  and  P,  or  fP,  are 
determined,  the  corresponding  stress  resultants  and  moments  are  readily 
calculated. 

We  now  present  two  classes  of  exact  solutions  in  three-dimensional 
isotropic  elasticity  theory.  We  consider  a  homogeneous  layer  of  uniform 
thickness  2h  with  its  mid-plane  at  Z=0.  It  is  convenient  to  employ  scaled 
variables  (x,  y,  z)  and  (u,  v,  w)  defined  as 

(x,  y,  z)  =  {X/a,  Y/a,  Zja),  (u,  v,  w)  =  (U/a,  Vja,  Wla'\  (2.6) 

where  a  is  a  characteristic  in-plane  dimension.  We  also  define  the  aspect 
ratio  E  as 

e=hla  (2.7) 

In  practice  it  is  expected  that  normally  e«  1,  but  the  analysis  is  exact  and 
does  not  depend  on  e  being  small. 

Solution  1 

Suppose  that  Mo(x,  y),  Vq{x,  y)  satisfy  the  classical  thin  plate  equations  for 
some  (unspecified)  value  of  7,  and 


^0  ~  ^O.x  ^O.y’  ^0  ~  ^0,x  ^O.y 


(2.8) 
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Thus 


2(1}  +  l)Ao.,  -fto., =0,  +  l)Ao.,  +  fio.x  =  0  (2-9) 

and  so 


V^Ao=0,  V^fio  =  0  (2.10) 

Then  the  three-dimensional  elasticity  equations  in  the  layer  with  elastic 
constant  tj  are  satisfied  by  u(x,  y,  z),  t!(x,  y,  z),  vv(x,  y,  z\  where 

o.r  I 

O.iJj 

(2.11) 

w(x,y,z)=  -Eri(z  +  Si)Ao 
where  5,,  Sj.  ^3  are  arbitrary  constants. 


■u(x,y,z)'  _  'uo(x,yy 
_t;(x,y,z)J  Lt>o(^.3’)J 


-(-S,z-hS3)<^(;;-h2) 


)M4-“ 

LAo.J  L 


(2.12) 


Solution  2 

Suppose  that  Wo(x,  y)  satisfies  the  biharmonic  equation 

?*wo  =  0 

Then  a  solution  of  the  three-dimensional  elasticity  equations  is 


(2.13) 

w(x,  y,  z)  =  WqIx,  y)  -f  £^>?(iz^  +  B^zA■B2)V^Wo 

where  Bi  to  B*  are  arbitrary  constants. 

For  either  solution,  the  stress  is  readily  calculated  by  substituting  (2.1 1) 
or  (2. 1 3)  into  the  stress-strain  relations. 

The  solutions  (2.11)  and  (2.13)  can  be  derived  in  various  ways;  one 
derivation  is  given  in  Ref.  1.  With  special  choices  of  the  constants  S,,  B,  they 
were  obtained  by  Michell  and  are  given  in  Love."*  However  these  authors 
do  not  seem  to  have  recognised  that  the  three-dimensional  solutions  are 
generated  by  two-dimensional  solutions  of  (2.9)  and  (2. 1 2). 

Let  us  now  consider  a  laminate  comprising  2N  -l- 1  layers  of  different 
isotropic  materials.  For  simplicity  we  consider  geometrically  symmetric 
laminates.  Quantities  relating  to  the  rth  layer  will  be  identified  by  the  index 
r;  the  layer  r = 0  contains  the  mid-plane  of  the  laminate  and  the  layer  r=N  is 
adjacent  to  the  upper  surface. 

Solutions  of  the  form  (2.1 1)  and  (2. 13)  are  adopted  in  each  layer,  with  the 
constants  Sj'’’  and  Bj'*  applying  in  the  rth  layer.  The  solutions  (2.1 1)  are 
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appropriate  for  stretching  deformations  and  (2.13)  for  bending  deforma¬ 
tions;  for  a  symmetric  laminate  these  solutions  uncouple  and  may  be 
treated  separately.  For  brevity  we  discuss  the  stretching  solution  only.  This 
has  available  the  3(^-1- 1)  constants  Conditions  which  have  to  be 
satisfied  are 

(a)  symmetry  conditions  at  the  mid-plane  2=0; 

(b)  continuity  of  displacement  and  traction  at  each  interlaminar  inter¬ 
face; 

(c)  traction-free  conditions  at  the  lateral  surfaces. 

It  is  found  that  all  of  these  conditions  can  be  satisfied  by  appropriate 
choices  of  Si'*,  Uo(x,  y)  and  y)-  Conditions  (a)  and  (b)  provide  a  set  of 
recurrence  relations  to  determine  the  Si'*;  except  for  a  single  normalising 
condition  these  constants  depend  only  on  the  laminate  stacking  geometry 
and  laminae  elastic  constants.  Hence  for  a  given  laminate  they  can  be 
calculated  once  and  for  all,  and  do  not  depend  on  any  particular 
boundary-value  problem  or  shape  of  plate.  The  displacements  Uq  and  Dq 
have  to  be  solutions  of  the  thin  plate  equations  (2.9)  with  the  elastic 
constant  tj  chosen  to  have  the  value  it  has  in  the  classical  laminate 
theory — that  is,  the  value  appropriate  for  the  equivalent  plate. 

We  may  therefore  follow  a  very  simple  procedure.  We  first  calculate  the 
elastic  constants  for  the  equivalent  plate,  as  in  classical  laminate  theory,  and 
the  laminate  constants  S)'*.  Then,  for  a  given  boundary  value  problem,  we 
solve  the  classical  laminate  theory  equations  to  determine  Uq  and  v^;  this 
may  be  done  by  any  available  method,  including  numerical  methods.  The 
displacement  in  each  layer  is  then  given  immediately  by  (2. 1 1 ),  and  the  stress 
by  substitution  in  the  stress-strain  relations.  The  bending  problem  may  be 
treated  similarly. 

One  of  the  quantities  of  interest  is  the  inter-laminar  shear  stress.  If  the 
traction  on  the  interface  between  layers  r—  1  and  r  is  then  we  find  that 

t,  =  -2^Hoeo(rio-^)-l-2  ^  grad  Aq  (2.14) 

where  fi„  rj,  relate  to  layer  s,  fj  relates  to  the  equivalent  plate,  e,  =  hja,  2h,  is 
the  thickness  of  the  layer  s,  and  Ag  arises  from  the  equivalent  plate  solution. 


3.  ANISOTROPIC  LAMINAE 

The  success  of  the  procedure  described  in  Section  2  arises  because  the 
right-hand  sides  of  (2. 1 1)  and  (2. 1 3),  regarded  as  power  series  in  z,  terminate 
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after  a  finite  number  of  terms,  and  reduce  to  polynomials  in  z.  When  the 
laminae  are  anisotropic,  this  situation  no  longer  obtains,  and  correspond¬ 
ing  closed  form  solutions  do  not  exist.  However  the  concept  of  generating 
three-dimensional  elasticity  solutions  from  two-dimensional  equivalent 
plate  solutions  remains  valid. 

We  consider  stretching  deformations  only.  In  classical  thin  plate  theory, 
the  governing  equations  for  the  mean  in-plane  displacement  components 
are 


Q 1 1  U.XX  +  2Qi6»^y  +  Q66  Uyy  +  Q16V.XX 
+  iQl2  +  Q66)^,xy  +  626  ^.yy  =  ^ 

(3.1) 

Q 1 6 “1" (6 1 2  “1"  0,16^, yy 

+  666  ^.xx  +  2626  •’.xji  +  Qzi^.yy  =  0 

where  Qi^  are  the  ‘reduced’  elastic  constants  of  the  material.  Solutions  of 
these  equations  can  be  represented  in  the  form  (e.g.  Leknitskii*) 

|^^j=2Re  {Pi4>i(x-(-Siy)-l-p2(?)2(x-l-S2y)}  (3.2) 

where  s,-  and  p,  are  the  eigenvalues  and  eigenvectors  of  the  eigenvalue 
problem  which  arises  by  seeking  solutions  of  (3. 1 )  of  the  form  u  =  p(p(x  +  sy). 

For  an  anisotropic  homogeneous  layer  we  seek  three-dimensional 
solutions  of  the  form 


u  =  U<jf){x-(-sy+et(r  — d)}  (3.3) 

where  u =(u,  v,  w)^.  On  substituting  (3.3)  into  the  stress-strain  relations  and 
the  equilibrium  equations,  and  provisionally  treating  s  as  known,  there 
results  an  eigenvalue  problem  with  eigenvalues  ±t,  and  corresponding 
eigenvectors  U,*  (i  =1,2, 3).  Thus  we  obtain  solutions  for  the  rth  layer  of  the 
form 

3 

u=  X  [Kr^U}^>>{x-l-sy-l-e,tr(z-dr^)} 

i=  1 

-I-  Kl'>-Ui'>-<p{x  -I-  sy-e,tl'\z  -  d)''- )}]  (3.4) 

where  K f’*  and  df’*  are  constants. 

For  a  laminate  we  assume  solutions  of  the  form  (3.4)  in  each  layer.  The 
constants  are  chosen  so  that  the  arguments  of  cp  are  continuous  at  each 
interface.  The  conditions  of  continuity  of  displacement  and  traction  at  the 
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interfaces  then  determine  to  leading  order  in  e  except  for  normalising 
factors.  It  remains  to  determine  the  parameter  s.  This  is  fixed  by  the 
traction-free  boundary  conditions  on  the  lateral  surfaces  of  the  plate.  It  is 
found  that  these  are  satisfied  to  leading  order  in  e  if  s  is  chosen  to  have  the 
value  given  by  classical  plate  theory,  for  the  equivalent  plate.  Thus,  we 
obtain  solutions  of  the  form 

u  =  Re  i  i  + 

(3.5) 

+  K  g'  -  UiJ’  -tpj{x  +  sjy-  £,tW(z  -  dr  - )}] 

where  Sj  are  determined  by  the  elastic  constants  for  the  equivalent  plate, 
(Pj{x  +  Sjy)  represent  a  solution  of  classical  laminate  theory  for  the 
equivalent  plate,  and  dj'’ *,  K  jj*  *,  *  can  be  determined  in  terms  of  the 
laminate  geometry  and  the  elastic  constants  of  the  laminae. 

Solutions  of  the  form  (3.5)  yield  exact  solutions  of  the  linear  elasticity 
equations  in  each  lamina.  However  they  satisfy  the  required  continuity 
conditions  at  the  inter-laminar  interfaces  and  the  traction-free  boundary 
conditions  at  the  lateral  surfaces  only  to  leading  order  in  £.  To  correct  this 
we  note  that  the  discrepancy  is  of  order  etp'jix  +  Sjy),  and  that  if  tpjlx  +  Sjy) 
generates  a  solution  of  the  governing  equations,  so  do  (pf’{x  +  Sjy),  where 
(pf^  denotes  the  nth  derivative.  Thus,  finally,  we  are  led  to  solutions  of  the 
form 

«  =  i  i  t  + 

(3.6) 

+  K  W  -  Ui5  -  Vf  {x  +  Sjy-  e,tW(z  -  d^  ' )}] 

Thus,  given  a  classical  thin  plate  solution  for  the  equivalent  plate,  the 
three-dimensional  solution  for  the  laminate  can  be  constructed  to  any 
required  order  in  e. 


4.  EDGE  EFFECTS 

The  solutions  outlined  in  Sections  2  and  3,  which  we  term  interior  solutions, 
are  exact  solutions  of  the  appropriate  elasticity  equations,  and  satisfy  the 
lateral  surface  boundary  conditions.  However  they  satisfy  edge  boundary 
conditions  only  in  an  average  sense.  To  satisfy  edge  boundary  conditions 
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point  wise,  it  is  necessary  to  superpose  additionaJ  solutions  which  neutralise 
the  mismatch  between  the  actual  boundary  conditions  and  the  boundary 
values  of  the  interior  solution.  It  is  to  be  expected  that  these  edge  solutions 
will  decay  exponentially  with  distance  from  the  edge;  however,  especially 
for  highly  anisotropic  materials  it  is  possible  that  the  exponential  decay  rate 
will  be  slow  and  the  edge  effects  may  penetrate  a  substantial  distance  into 
the  plate.  Also,  the  edge  solutions  may  include  stress  singularities  at  the 
intersection  of  an  edge  with  an  interlaminar  interface. 

Although  the  determination  of  the  edge  solution  is  a  separate  problem, 
the  interior  solution  can  give  useful  information  about  the  magnitude  of  the 
stress  in  the  edge  solution.  In  particular,  the  magnitude  of  the  mismatch 
between  the  actual  boundary  conditions  and  the  boundary  values  of  the 
outer  solution  controls  the  magnitude  of  the  edge  solution.  It  is  tentatively 
suggested  that  some  suitable  measure  of  this  mismatch  might  serve  as 
a  design  parameter  which  could  be  used  as  a  measure  of  susceptibility  to 
edge  delamination. 
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ABSTRACT 

We  review  some  of  the  recent  mathematical  progress  on  the  effective  moduli  of 
composites.  Specific  attention  is  devoted  to  mathematically  precise  definitions 
of  effective  moduli,  new  methods  for  bounding  effective  moduli,  new  construc¬ 
tions  of  mixtures  with  explicitly  computable  properties,  and  applications  to 
structural  optimization. 


1.  INTRODUCTION 

We  are  concerned  with  materials  that  are  spatially  heterogeneous  on 
a  suitably  small  length  scale,  and  with  linear  models  of  material  behavior, 
for  example  linear  elasticity.  The  effective  moduli  of  such  a  ‘composite’ 
describe  its  overall,  large-scale  behavior.  They  have  long  been  an  object  of 
study  by  physicists  and  materials  scientists;  selective  reviews  of  the 
extensive  literature  include  Refs.  14, 22, 68, 69, 72.  More  recently,  the  study 
of  effective  moduli  has  attracted  the  attention  of  a  growing  community  of 
mathematicians  as  well.  The  relatively  new  notions  of  homogenization  and 
G-convergence  provide  a  firm  mathematical  foundation;*®  ®^'®*-’*  more¬ 
over,  the  effective  moduli  of  composites  have  been  linked  to  fundamental 
issues  arising  in  the  optimal  control  of  certain  distributed  parameter 
systems,  and  to  deep  questions  involving  the  lower  semicontinuity  of 
variational  functions,  see,  e.g.  Refs.  1, 12, 28, 30, 32, 38, 39, 49, 50, 58, 67.  The 
specific  questions  about  effective  moduli  raised  by  these  new  applications 
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are  sometimes  difTerent  from  those  that  were  the  focus  of  the  older 
literature:  for  example,  applications  to  structural  optimization  require  the 
specification  of  all  (anisotropic)  composites  attainable  as  mixtures  of  given 
components  in  specified  proportions.  However,  the  mathematical  tools 
developed  to  address  such  questions  have  also  led  to  new  results  that  are 
very  much  within  the  purview  of  the  older  theory.  Examples  include  the 
simultaneous  attainability  of  the  Hashin-Shtrikman  shear  modulus  and 
bulk  modulus  bounds;'^'^*’'^^**  the  validity  of  a  conjecture  of  Schulgasser 
about  the  effective  conductivity  of  polycrystalline  composites;*  and  the 
attainability  of  certain  mean  field  theories.^-^* 

The  goal  of  this  paper  is  to  review  selected  aspects  of  this  recent 
mathematical  progress,  which  it  is  hoped  will  be  of  interest  to  a  bro  .d 
community  of  specialists  in  materials  science.  It  should  be  emphasized  that 
the  ideas  presented  here  are  a  synthesis  of  the  work  of  many  individuals,  and 
that  the  selection  of  topics  is  somewhat  arbitrary — in  no  way  representing 
a  comprehensive  survey  of  the  most  important  recent  developments. 


2.  MATHEMATICALLY  PRECISE  DEFINITIONS 
OF  EFFECTIVE  MODULI 

We  are  concerned  with  mixtures  of  continua  on  a  length  scale  small 
compared  to  that  on  which  the  loads  and  boundary  conditions  vary,  but 
still  large  enough  for  continuum  theory  to  apply.  Such  a  ‘composite’  is 
clearly  an  idealization:  it  represents  the  limiting  behavior  of  a  sequence  of 
structures,  as  the  ratio  e  =  IjL  relating  the  ‘microscopic’  length  scale  /  to  the 
‘macroscopic’  one  L  tends  to  zero.  There  are  in  fact  several  distinct  theories, 
differing  as  to  the  form  assumed  for  the  fine  scale  structure.  A  periodic 
composite  is  one  whose  microscopic  structure  is  periodic  with  a  specified 
unit  cell;  a  random  composite  is  one  whose  fine  scale  structure  is  a  stochastic 
process  with  specified  statistics.  There  is  also  a  third  approach  which  makes 
no  such  hypothesis  on  the  fine  scale  structure,  appealing  instead  to 
a  compactness  theorem  for  systems  of  partial  differential  equations.  This 
last  theory,  known  variously  as  G-convergence  or  homogenization,  re¬ 
presents  in  a  sense  the  most  general  approach. 

To  fix  ideas,  let  us  focus  the  discussion  on  mixtures  of  two  isotropic, 
linearly  elastic  materials  in  R‘'(d = 2  and  d=3  being,  of  course,  the  cases  of 
physical  interest).  Each  of-the  component  materials  is  characterized  by 
a  bulk  modulus  k,  and  a  shear  modulus  (i=  1,2),  determining  a  unique 
Hooke’s  law  tensor  A, — ^a  symmetric  linear  map  on  the  space  of  symmetric 
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Aie=Kt(tTe)I  +  2fiJie-  j(tr  f!)I 
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(2.1) 


for  any  symmetric  tensor  e.  The  associated  ‘elastic  energy’  quadratic  form  is 
the  inner  product  of  stress  and  strain: 


(A/e,  e)= 


(tr  + 


(2.2) 


A  structure  which  mixes  the  two  materials  will  have  a  spatially  varying 
Hooke’s  law,  equal  to  either  A  i  or  ^2  at  each  material  point  x.  Introducing 
a  parameter  e,  representing  (at  least  in  the  periodic  and  random  cases)  the 
length  scale  of  the  microstructure,  the  spatially  varying  Hooke’s  law  is 


A%x)=xMx)Ai+xMx)A2  (2.3) 

where 


Zf(^)  = 


on  the  set  occupied  by  material  i 
elsewhere 


(2.4) 


so  that  X2  =  ^—X\-  By  definition  the  structure  is  periodic  (with  cubic) 
symmetry)  if 

Z‘(Jc)  =  ):i(i)  for  some  function  Xi(y), 

taking  only  the  values  0  and  1,  defined  for  all  y  e  R** 

and  periodic  in  each  component  of  y  with  period  1.  (2.5) 

An  example  would  be  a  periodic  array  of  spherical  inclusions  centered  on 
a  cubic  lattice  of  mesh  e,  each  sphere  having  radius  ep{p<i).  In  the  random 
case  there  is  an  additional  variable  cu,  belonging  to  a  suitable  probability 
space: 

Xf (x,  (u)  =  Xi(f ,  (o)  for  some  stochastic  process  a>->x,(  y,  cu), 
defined  for  y  e  R"*  and  « in  a  probability  space,  and  taking 
only  the  values  0  and  1.  It  is  required  that  x,  be  translation 
invariant,  in  the  sense  that  tu-^Xify  +  c,  w)  gives  the 
same  stochastic  process  for  each  c  e  R"*.  Furthermore,  the 
translations  are  assumed  to  be  ergodic,  so  that  ensemble 
averaging  is  equivalent  to  spatial  averaging.  (2.6) 

An  example  would  be  a  family  of  (possibly  overlapping)  spherical  inclusions 
of  radius  sp  whose  centers  have  a  multidimensional  Poisson  distribution. 
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the  expected  number  of  balls  in  a  unit-sized  region  being  of  order  The 
hypotheses  (2.5)  or  (2.6)  specify  rather  precisely  the  character  of  the  fine 
scale  structure.  The  G-convergence  approach,  by  contrast,  makes  no  such 
hypothesis: 

/ifx)  is  any  family  of  functions  taking  only  the  values  0 
and  1,  parametrized  by  £-»0,  and  iifx  =  1  —x\  •  (2.7) 

It  is  specifically  not  assumed  in  (2.7)  that  e  represents  the  length  scale  of  the 
microstructure:  even  a  sequence  which  has  no  clear  separation  of  scales  is 
permitted.  Clearly  (2.7)  includes  both  the  periodic  case  and  the  random  one; 
indeed,  in  our  opinion  it  includes  any  reasonable  notion  of  a  linearly  elastic 
composite  obtained  by  mixing  two  materials  (with  perfect  bonding  at  all 
material  interfaces). 

The  tensor  of  effective  moduli  A*  is  simply  the  Hooke’s  law  tensor  of  the 
composite.  It  represents  the  limiting  behavior  of  the  mixture  as  £->0.  This 
means  that  for  any  (£-independent)  load  f  the  associated  elastostatic 
displacement  — which  solves  the  equilibrium  equations 

<7*  = /IV 


div  (f—f 

with  an  appropriate  boundary  condition — converges  as  £-+0  to  «*,  the 
solution  of  the  corresponding  system  with  A‘  replaced  by  A*.  The  starting 
point  of  the  mathematical  theory  is  the  existence  of  effective  moduli.  In  the 
spatially  periodic  and  stationary  stochastic  contexts  (2.5),  (2.6),  translation 
invariance  assures  that  the  tensor  .4*  of  effective  moduli  is  constant.  For 
periodic  composites  it  can  be  given  in  terms  of  the  solutions  of  certain 
canonical  ‘cell  problems’,  see  for  example.  Refs.  8,  60,  but  we  prefer  this 
variational  characterization,  cf.  Ref.  64: 

(/l*^a=inff  {MyM  +  e{<f>)l  ^ +  €{(!,)) dy  (2.9) 

*  JQ 

in  which 

My>Xiiy)'^i+X2(y)'^2  (2io) 

Q  =  [0*  1]“*  the  unit  cell  of  the  perodic  structure,  (p  varies  over  Q-periodic 
displacement  fields,  and  e{<l>)=]{V(f>  +  V<p^)  is  the  linearized  strain  as- 
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sociated  to  tj).  An  entirely  analogous  formula  is  available  in  the  random 
case,  cf.  Refs.  19,  33,  55,  70: 

inf  E[(^e,e)]  (2.11) 

in  which  E  represents  the  ensemble  average  and  e  ranges  over  stationary, 
random  strain  fields  with  mean  value  In  the  more  general  G-convergence 
setting  (2.7)  there  is  no  hypothesis  of  translation  invariance,  so  the  tensor  of 
effective  moduli  A*{x)  can  vary  with  x.  Moreover,  there  is  obviously  not 
enough  structure  to  give  a  formula  as  explicit  as  (2.9)  or  (2.11).  But  it  is 
nevertheless  true  that  for  any  sequence  xf  as  in  (2.7)  there  is  a  subsequence 
£'-♦0  for  which  there  exists  a  limiting  tensor  of  effective  moduli  /l*(x),  see 
for  example  Refs.  48,  62,  65,  71. 

We  shall  be  interested  in  bounds  for  A*  in  terms  of  the  volume  fractions 
of  the  component  materials,  so  let  us  note  here  how  to  express  these  volume 
fractions  in  each  of  the  different  settings.  For  the  periodic  composite  (2.5) 
the  volume  fraction  of  material  i  is  the  proportion  of  the  period  cell 
occupied  by  it: 

0.=  f  (2.12) 

Jo 

Similarly,  in  the  stationary,  random  case  (2.6)  it  is  the  expected  value  of 

0.  =  E(z,)  (2.13) 

In  the  G-convergence  context  (2.7)  it  is  instead  given  by  the  L® — weak*  limit 

0,(x)=wk*limxf(^)  (2.14) 

no  longer  necessarily  constant,  defined  by  the  property  that 

|xi(^)a(.>c)dx-»j0i(x)g(x)dx  (2.15) 

for  continuous  functions  g. 

These  notions  of  effective  moduli  are  easily  seen  to  be  equivalent  to  the 
operational  definitions  more  commonly  used  in  materials  science,  based  on 
the  average  stress  and  strain  or  average  elastic  energy  in  a  physical  domain 
that  is  large  compared  with  the  microstructure  but  small  compared  with  the 
length  scale  of  the  loads  and  boundary  conditions,  sec  for  example.  Refs.  22, 
24.  They  are  important  for  the  development  of  a  proper  mathematical 
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theory,  because  they  make  it  possible  to  give  fully  rigorous  proofs  of  results 
about  effective  moduli.  But  why  should  they  be  of  interest  to  a  materials 
scientist?  One  answer  lies  in  the  following  ‘density’  result;'*  if  an  algebraic 
relation  between  the  tensor  of  effective  moduli  and  the  component  volume 
fractions  holds  for  all  spatially  periodic  composites  (or  for  all  stationary, 
stochastic  composites),  then  it  holds  in  the  more  general  context  of 
G-convergence  as  well.  Thus,  for  bounds  on  effective  moduli  in  terms  of 
volume fractions  alone,  neither  long-range  disorder  nor  a  definite  separation  of 
scales  is  relevant.  This  resolves  a  point  which  has  been  the  object  of 
considerable  controversy  in  the  literature,  see  for  example,  Ref  22. 


3.  NEW  METHODS  FOR  BOUNDING  EFFECTIVE  MODULI 

A  typical  goal  of  the  new  mathematical  theory  is  the  so-called  G-closure 
problem:  find  the  precise  set  of  Hooke’s  laws  A*  achievable  by  mixing  two 
given  isotropic,  elastic  materials  in  specified  proportions.  The  motivation 
comes  from  applications  to  structural  optimization,  as  we  shall  explain  in 
Section  5.  The  special  case  when  is  isotropic  was  considered  by  Hashin 
and  Shtrikman,^^  under  the  further  hypothesis  that  the  component  mater¬ 
ials  are  well-ordered  i.e.  that 

(3.1) 

They  gave  upper  and  lower  bounds  for  the  effective  bulk  and  shear  moduli, 
K*  and  p*,  which  are  now  known  to  be  simultaneously  achievable. An 
improvement  of  the  Hashin-Shtrikman  bounds  can  be  found  in  Refs.  10 
and  47,  but  the  precise  set  of  attainable  isotropic  composites  is  still  not 
known.  In  any  event,  results  of  this  kind— concerning  A*  with  specified 
symmetry — are  not  adequate  for  applications  to  structural  optimization, 
since  the  best  composites  for  use  in  an  optimal  structure  can  (and  generally 
will)  be  fully  anisotropic.  While  the  complete  solution  of  the  G-closure 
problem  seems  beyond  the  reach  of  current  methods,  the  analogues  of  the 
Hashin-Shtrikman  bounds  on  k*  and  p*  are  now  understood  for  fully 
anisotropic  composites.^  '*  *’  In  particular,  we  now  know  those  parts  of  the 
boundary  of  the  G-closure  which  represent  the  ‘strongest’  and  the  ‘weakest’ 
anisotropic  composites. 

In  the  course  of  exploring  these  and  other  bounds  for  effective  moduli, 
a  number  of  powerful  new  tools  have  been  introduced.  The  well-known 
Hashin-Shtrikman  variational  principles  have  been  applied  in  new 
ways,’  *'^*  ’*  *’  and  new  variational  principles  have  been  introduced. 
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obtained  from  more  classical  ones  by  the  addition  of  a  quadratic 
null-Lagrangian.*-^^  In  addition,  entirely  new  approaches  have  been 
introduced:  one  is  based  on  an  equivalence  between  bounds  for  effective 
moduli  and  the  lower  semicontinuity  of  certain  variational  func- 
tionals;^®'^^  ®^  another  uses  the  fact  that  the  effective  moduli  depend 
analytically  on  the  component  properties, a  third  uses  ‘com¬ 
pensated  compactness’  to  construct  certain  lower  semicontinuous  func¬ 
tionals,*  and  a  fourth  makes  use  of  Hilbert  space  decomposi¬ 

tions  and  continued  fractions.^^  (These  references  represent  a  mere 
sampling  of  the  relevant  literature  in  each  area.)  The  interested  reader  will 
find  several  of  these  new  methods  applied  to  a  single  problem  in 
a  self-contained  manner  in  Ref  27.  The  power  and  limitations  of  these 
various  methods  are  just  beginning  to  be  understood,  as  are  the  relation¬ 
ships  among  them.'*® 

To  convey  some  of  the  flavor  of  these  new  developments,  we  present  in 
detail  one  of  the  recently  established  bounds,  an  upper  bound  on  the  elastic 
energy  quadratic  form.  There  is  of  course  a  well-known  bound  due  to 
Paul:®® 


{A*U)^et{A,U)  +  d,{A,U)  (3.2) 

where  0,  is  the  volume  fraction  of  the  ith  material,  i  =  1,2.  This  bound  is 
sharp,  in  the  sense  that  for  certain  choices  of  the  ‘average  strain’  ^  there  is 
a  microstructure  whose  associated  A*  achieves  equality  in  (3.2).  However, 
for  most  choices  of  (3.2)  is  not  saturated  by  any  composite;  therefore 
a  better  bound 


(3.3) 

is  possible.  We  shall  in  fact  prove  the  optimal  bound  of  this  type,  in  other 
words  one  which  is  saturated,  for  each  by  an  appropriately  chosen 
mixture  of  the  two  given  materials.  The  method,  which  is  based  on  the 
Hashin-Shtrikman  variational  principle,  requires  that  the  component 
materials  be  well  ordered.  Our  presentation  follows  that  of  Ref  26; 
equivalent  results  can  be  found  presented  somewhat  differently  in  Refs.  3, 4 
and  45.  The  function  F  on  the  right  of  (3.3)  is  given  by  (3.16)  below,  as  the 
extremal  value  of  a  finite-dimensional  mathematical  programming 
problem. 

As  discussed  in  Section  2,  it  is  sufficient  to  prove  the  bound  for  spatially 
periodic  composites.  We  may  therefore  fix  Q  =  [0, 1]*  as  the  period  cell;  the 
microstructure  is  determined  by  the  indicator  functions  /ify)  and  X2(>')  = 
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1  —  Xiiy),  y^Q,  constrained  by  the  given  volume  fractions  (2.12);  and  the 
effective  Hooke’s  law  is  determined  by  (2.9). 

The  first  step  is  to  derive  the  Hashin-Shtrikman  variational  principle: 


^  +  e{(j)))Xidy 


+  (*  ((^2-'4,)  ‘ff,(T)xidy 
jQ 


+  \  (A^(i-ye((i>))A  +  emdy  (3.4) 

Jo 

for  any  ^-periodic  displacement  field  </>,  and  any  Q-periodic  field  of 
symmetric  tensors  a.  The  proof  is  elementary:  expanding  the  pointwise 
inequality 

+  (3.5) 

and  multiplying  by  Xi  gives 

-  2(0,  ^  +  e(<l>))xi+Xi(('42-Ai)~^(7,(7)  (3.6) 

The  left  side  equals 

{(A-A,)(i  +  e(<P)),^  +  e{c/>))  (3.7) 

therefore  integrating  over  Q  and  applying  (2.9)  we  conclude  (3.4). 

The  next  step  is  to  specialize  to  constant  a,  and  to  evaluate  the  integrals  in 
(3.4)  wherever  possible.  This  gives 

({A*-A2K,^)  +  20,(ff,O-0t((A2-A,r'<7,(7) 


* 

Jo 


e((^))dy+  (A2e{<p),e{<P))dy 

JQ 


for  any  Q-periodic  displacement  field  (f>. 

The  third  step  is  to  minimize  the  expression  on  the  right  over  <f>.  This 
amounts  to  solving  the  elastostatic  equilibrium  equation 

div(/l2^0))-div(<T;f,)  =  O  (3.9) 

with  a  periodic  boundary  condition.  It  is  convenient  to  use  Fourier  analysis: 
since  A  2  and  o  are  constant,  (3.9)  determines  the  Fourier  transform  of  4>  at 
each  frequency  keZ‘‘  directly  in  terms  of  the  transform  of  Xi  at  the  same 
frequency.  After  some  algebra,  one  finds  that  the  extremal  value  of  the  right 
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side  of  (3.8)  is 


- 1 

(3.10) 

where 

Xi(y)=  S 

(3.11) 

keZ-‘ 


and  for  any  unit  vector  v  e  R‘',y(f)  is  the  ‘degenerate  Hooke’s  law’  defined  by 


f{v)<T= 


dK2  +  2(d—\)fl2 
1 


{(TV,  v)v  o  V 


H - [(oi?)  o  v—((TV,v)v  o  v] 


(3.12) 


Here  <t  is  any  symmetric  tensor,  and  we  use  the  notation  vQw  = 
i(D®w  +  w(8)t))  for  the  symmetric  tensor  product  of  two  vectors  in  R"*. 

It  remains  to  eliminate  the  explicit  dependence  of  the  bound  on j ,  which 
is  after  all  arbitrary  except  for  the  volume  fraction  constraint.  We  use  this 
constraint  to  see  that 


(Xi-0i)'dy  =  0,0,  (3.13) 

JQ 

whence  by  Plancherel’s  theorem 

=  (3.14) 

This  gives  a  bound  on  the  ‘nonlocal’  term: 

(3.10)  <  —  0, 02  "tin  {f{v)<T,(T)  (3.15) 

101=1 

Substitution  into  (3.8)  gives  a  bound  on  A*  which  still  depends  on  the 
choice  of  a  symmetric  tensor  a,  and  minimization  over  a  gives  a  result  of  the 
desired  form  (>1*^,(^XF,  with 

F  =  (/I2  ^, ^)  +  0,  •  min  { - 2(<t, ^) + ((542  - .4 , ) -  >  ff, (7) 

a 

-02  min  {f{v)<T,(T)}  (3.16) 

|o|  =  l 

Our  interest  in  this  bound  lies  in  the  fact  that  it  is  the  best  possible  bound 
for  (A*  (J,  (J)  in  terms  of  the  given  parameters  0, ,  02  =  1  —  0, ,  and  the  bulk 
and  shear  moduli  of  the  component  materials  /c,  <k:2,/i,  </i2-  This  will  be 
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proved  in  the  next  section,  as  an  application  of  the  formula  for  the  effective 
behavior  of  a  sequentially  laminated  composite. 


4.  CONSTRUCTION  OF  MIXTURES  WITH  EXPLICITLY 
COMPUTABLE  EFFECTIVE  MODULI 

For  most  microstructures  there  is  no  explicit,  algebraic  formula  for  the 
tensor  of  effective  moduli  A*;  one  must  make  do  instead  with  a  variational 
principle  such  as  (2.9)  or  (2.11),  or  perhaps  with  the  partial  differential 
equation  characterizing  its  extremal.  If  this  were  the  only  available  tool  it 
would  be  virtually  impossible  to  establish  the  optimality  of  any  bound! 
Fortunately  there  are  certain,  rather  special  microstructures  for  which  the 
effective  moduli  are  computable;  and,  remarkably,  this  class  of  composites 
is  rich  enough  to  demonstrate  the  optimality  of  a  variety  of  bounds, 
including  (3.3). 

There  are  some  simple  and  more  or  less  classical  examples  of  composites 
with  explicitly  computable  properties.  One  example  is  that  of  a  layered 
microstructure;*' another  is  the  ‘concentric  sphere  construction’, 
which  was  used  by  Hashin  in  Ref.  73  to  prove  the  optimality  of  their  bulk 
modulus  bounds.  It  is  natural  enough  to  iterate  such  constuctions,  for 
example  layering  together  two  composites  each  of  which  has  its  own 
fine-scale  structure,  obtained  perhaps  by  layering  or  by  a  version  of  the 
concentric  sphere  construction.  This  idea,  which  can  be  found  in 
Bruggeman’s  work,“  has  been  rediscovered  by  various  individuals  and 
applied  to  prove  the  attainability  of  many  different  bounds,  e.g.  Refs.  3-5, 
17,  18,  26,  34,  35,  37,  38,  42,  61,  66. 

An  important  new  development  concerns  the  attainability  of  certain 
mean  field  theories.  The  formulas  they  predict  for  the  tensor  of  effective 
moduli  A*  were  originally  intended  as  approximations,  not  as  exact  results. 
Nevertheless,  it  has  recently  been  shown  that  certain  effective  medium 
theories  are  exactly  attainable  by  composites  with  appropriately  chosen 
microstructures.^'^*'**  Obviously,  this  result  greatly  expands  the  class  of 
composites  with  explicitly  computable  effective  moduli — particularly  since 
these  effective  medium  theories  (the  ‘coherent  potential  approximation’  and 
the  ‘differential  effective  medium  theory’)  have  been  widely  studied  in  the 
mechanics  literature,  see  e.g.  Refs.  74,  75. 

The  microstructures  that  arise  from  these  constructions  are,  it  should  be 
understood,  somewhat  idealized  materials.  They  are  higly  ordered,  neither 
periodic  nor  stochastic  in  character,  and  they  frequently  involve  multiple 
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length  scales.  It  may  seem  like  cheating  that  we  allow  the  use  of  such 
microstructures  to  establish  the  attainability  of  a  bound,  whereas  the  proof 
of  the  bound  may  make  use  of  a  special  structure  such  as  periodicity.  This  is 
in  fact  perfectly  legitimate;  indeed,  it  is  here  that  we  use  the  power  of  the 
mathematical  theory.  The  point  is  that  these  constructions  fit  perfectly  into 
the  mathematical  context  of  G-convergence  (see  especially  Ref.  2);  therefore, 
by  the  ‘density’  result  mentioned  at  the  end  of  Section  2,  their  effective 
moduli  can  be  approximated  arbitrarily  well  by  those  of  spatially  periodic 
composites.  Actually,  it  is  quite  natural  to  use  the  most  restrictive  possible 
setting  for  proving  bounds,  and  the  most  general  one  for  showing  that  they 
are  achieved. 

The  remainder  of  this  section  is  devoted  to  a  discussion  of  sequentially 
laminated  composites,  and  to  a  proof  of  the  attainability  of  the  new  upper 
bound  (3.3).  Closely  related  ideas  and  results  can  be  found  in  Refs.  3, 4, 26, 
45.  The  construction  of  a  sequentially  laminated  composite  is  an  iterative 
procedure,  producing  a  microstructure  that  has  several  different  length 
scales.  A  laminar  composite  of  rank  1  is  obtained  by  layering  two  initially 
given  materials,  specifying  the  proportion  of  each  and  the  layer  direction, 
and  using  a  small  parameter  e,  as  the  layer  thickness.  As  e,  ->0,  the  elastic 
behavior  is  described  by  an  effective  Hooke’s  law  Cj .  A  laminar  composite  of 
rank  2  is  obtained  by  layering  two  laminar  composites  of  rank  1,  again 
specifying  the  proportion  of  each  and  the  layer  direction,  and  using  another 
small  parameter  Sj  for  the  layer  thickness.  As  e,,e2->0  with  ei«e2,  the 
elastic  behavior  is  described  by  an  effective  Hooke’s  law  C2 .  This  process 
can  clearly  be  continued  indefinitely:  the  general  sequentially  laminated 
composite  of  rank  r  is  obtained  by  layering  together  two  sequentially 
laminated  composites  of  rank  r-  1.  We  shall  consider  here  only  a  special 
case,  in  which  one  of  these  two  materials  is  the  isotropic  one  with  shear 
modulus  P2  modulus  K2  at  each  successive  stage.  An  elegant,  iterative 

formula  for  representing  the  effective  moduli  of  such  a  composite  was 
derived  in  Ref.  1 7,  following  a  method  developed  for  scalar  equations  in  Ref. 
66.  We  now  give  a  derivation  of  this  result. 

The  basic  building  block  is  a  formula  for  the  effective  tensor  C  cor¬ 
responding  to  a  layered  mixture  of  the  isotropic  material  with  Hooke’s  law 
A  2  and  a  general  elastic  material  with  Hooke’s  law  B,  using  layers 
orthogonal  to  the  unit  vector  ceR",  and  using  proportions  P2 
pg=  I— P2  of  ^2  resp)ectively: 

Pg(A2-C)-'tT=(A2-Br'o-p2f(v)o  (4.1) 

for  any  symmetric  tensor  <7.  Here  f(v)  is  the  same  degenerate  Hooke’s  law 
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that  arose  in  our  proof  of  the  bound,  defined  by  (3.12).  In  writing  (4.1)  we 
have  implicitly  assumed  that  A-^  —  C  and  /Ij — are  invertible,  when  viewed 
as  symmetric  linear  maps  on  the  space  of  symmetric  tensors.  This  is  the  case 
whenever  since  then  C<A2  as  well,  by  Paul’s  bound  (3.2);  this 

ordering  hypothesis  will  be  sufficient  for  our  purposes,  since  we  are 
concerned  with  mixtures  of  two  well-ordered  isotropic  materials,  i.e.  (3.1) 
holds.  (There  is  a  version  of  (4.1)  without  invertibility  hypotheses,  see  for 
example  Ref.  17.)  To  prove  (4.1),  one  must  of  course  begin  with  a  charac¬ 
terization  of  C.  In  a  layered  composite  of  the  type  under  consideration,  the 
local  values  of  the  stress  and  strain  are  essentially  constant  within  each 
component.  Therefore,  arguing  for  example  as  in  Ref.  40,  the  calculation  of 
given  ^  is  easily  reduced  to  this  algebraic  problem:  find  a  pair  of 
symmetric  matrices  (representing  the  strain  in  the  layers  occupied 

by  materials  A  2  and  B  respectively)  such  that 

P2  Pb^b~  ^ 

ijg  — (J2  =  y  O  w  for  some  vveR"  (4.2a-c) 

(A2^2-BiB>  =  ^ 

The  first  relation  says  that  is  the  average  strain;  the  second  is  the 
consistency  condition  for  the  existence  of  a  deformation  with  the  specified 
piecewise  constant  strain  (recall  that  t?0>v=(j><8>  w-i-w®  t>)/2);  and  the 
third  represents  the  continuity  of  the  norma)  stress  at  the  layer  interface.  In 
terms  of  these  quantities,  is  determined  by 

a  =  -42  ^2  + Pi.  5^1.  (4-2d) 

which  identifies  it  as  the  average  stress.  The  solution  of  (4.2a-d)  is  easiest 
to  represent  in  terms  of  a=(A2  —  C)^.  One  calculates  that  and  are 
given  in  terms  of  a  by 

^B  =  PB\Ai-B)'^a,  =  (4-3) 

where  w  €  R"  is  chosen  so  that 

P«  ^4  2(1;  O  w)  =  2(ai.)  O  i. — (<Ti;,  t.)t!  ©  1;  (4.4) 

whence 

Pi»['42(‘’  O  w)]d=<ti;  (4.5) 

The  unique  w  satisfying  (4.4)  is 

=  Pfl  M  1 — ,  v)v  +  ~(av  -  {(TV,  v)v) 

ldK2-k-2id-l)ti2  P2 


(4.6) 
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and  it  has  the  property  that 

PBvQw=f(v)a  (4.7) 

with  f(v)  defined  by  (3.12).  Therefore 

=  Pb‘(^2-B)'‘<7-Pb (4-8) 

which  is  precisely  the  desired  formula  (4.1). 

Now  consider  a  sequence  Q.Ci.Cj,...  of  effective  tensors  such  that 

Co  =  Ai  represents  an  isotropic  material  with  bulk 

modulus  fCj  and  shear  modulus  fii,  (4.9a) 

and,  for  1, 

C,  is  obtained  by  layering  A  2  with  C,_,  in  volume 
fractions  a,  and  (1  -  a,)  respectively,  using  the  unit 
vector  V,  as  the  layer  normal.  (4.9b) 

Evidently,  C,  represents  the  effective  behavior  of  a  certain  sequentially 
laminated  composite  of  rank  r.  The  volume  fraction  of  .4^  in  C,  is 

'•^1;  Po  =  0  (4.10) 

1 

A  formula  for  C,  is  easily  obtained  by  iterating  (4.1): 

(1-^,)(A,-Q-‘=(A2-A,)-‘-  t  (4.11) 

i=  1 

Let  us  terminate  this  process  at  r  =  yv,  and  write 
0^  =  =  overall  volume  fraction  of 

0^  =  \—Pf^  =  overall  volume  fraction  of  A^  (4.1 2) 

A  *  =  C^= effective  Hooke’s  law  of  the  associated  rank  N 
composite. 

It  is  easy  to  see  that  the  sequence 

m..=(^,-^,_,)/)?^,  (4.13) 

can  be  any  nonnegative  sequence  which  sums  to  1,  by  making  an 
appropriate  choice  of  the  parameters  {a,}.  Thus  we  have  shown  that  for  any 
integer  N^l,  any  unit  vectors  in  R"*,  any  real  numbers  = ,  with 

0<m,^  1  and  Snj,  =  1  and  any  real  number  02,O<02  <  E  there  is  a  sequen- 
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tially  laminated  composite  made  by  mixing  A  j  and  A  2  as  in  (4.9),  using  overall 
volume  fractions  0^  =  1— $2  and  $2  respectively,  whose  effective  Hooke’s  law 
A*  is  characterized  by 

e,(A2-A*)-^=(A2-A,r^-e2  j^mM  (4.14) 

i=i 

We  now  apply  this  construction  to  establish  the  optimality  of  the  new 
upper  bound  (3.16)  Our  task  is  to  show  that  for  each  symmetric  tensor 
there  is  a  choice  of  the  parameters  such  that  A*,  defined  by  (4.14), 

satisfies  (A*^,<^)=F  with  F  as  in  (3.16).  Now,  (3.16)  gives  F  in  terms  of 
a  mathematical  programming  problem 

min{-2(<T,^)+((.42-y4,)”‘(T,(T)-02  min  (/(t))ff,<T)}  (4.15) 

"  i«l  =  i 

over  symmetric  tensors  a,  so  it  is  reasonable  to  expect  the  proper  choices  of 
to  emerge  from  the  optimality  conditions  for  (4.15).  Since  the  last 
term  is  not  a  smooth  function  of  o,  it  is  natural  to  use  the  methods  of 
‘nonsmooth  analysis’,  see  for  example  Ref  1 5.  To  this  end  we  rewrite  (4. 1 5) 
as 

min  { -  2(<t,  +  g(a) }  (4.16) 

with 

g(<7)  =  max  ((^  2  “  -^ , ) '  *  -  02f(v)<r,  ff)  (4. 1 7) 

For  each  fixed  v  the  expression  on  the  right  is  a  positive,  quadratic  function 
of  a  (one  way  to  establish  positivity  is  to  make  use  of  (4.1)).  Therefore  g  is 
convex,  and  the  optimality  condition  for  (4.17)  is  that  for  any  extremal  or* 

2ieeg((T*)  (4.18) 

where  dg{a*)  is  the  subdifferential  of  g  at  a*  (see  for  example  Ref  15, 
2.3. 1-2.3.3  and  Corollary  1,  §2.3).  Moreover,  dg{(T*)  is  the  closed  convex  hull 
of  the  differentials  of  the  various  quadratic  forms  in  (4. 1 7)  as  r  ranges  over 
all  extremals  (see  for  example.  Ref  15,  §2.8,  Corollary  1).  Since  the  space  of 
symmetric  tensors  is  finite  dimensional,  each  element  of  the  closed  convex 
hull  is  in  fact  a  convex  combination  of  finitely  many  extreme  points. 
Therefore  the  optimality  condition  (4.18)  becomes 

^  =  (/l2-/l,)-‘<r*-02  Z  m.ffvM* 

i=  1 


(4.19) 
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with  Im—l,  |t),l  =  l,  N<oo,  and 

g{a*)=({A2-Ai)~^a*,(T*)-d2(f{Vi)a*,a*),  (4.20) 

Comparing  (4.19)  with  (4.14),  we  see  that 

^  =  (4.21) 

where  A*  is  the  sequentially  laminated  composite  of  rank  N  constructed 
using  (ntitVijiLi.  We  claim  that  this  A*  satisfies  (A*^,^)=f.  Indeed,  the 
value  of  F  is 

F  =  (A2U)+0A-  2(«t*,  a  +  (4.22) 

using  (3.16)  and  the  fact  that  <r*  k  extremal  for  (4.15).  We  have 

(a*,0  =  g(tx*)  (4.23) 

by  (4.19)  and  (4.20),  so  (4.22)  becomes 

F=(A2U)-0i(<r*,^)  (4.24) 

But  01  ff*  =  (A2  —  A*)^  by  (4.21),  and  substitution  gives  the  desired  result 

F=(A*^,a 


5.  APPLICATIONS  TO  STRUCTURAL  OPTIMIZATION 

The  recent  interest  in  optimal  bounds  on  the  effective  moduli  of  composites 
has  been  stimulated  in  large  part  by  applications  to  structural  optimization, 
see  for  example.  Refs.  1, 28,  38,  39, 49,  50,  67.  That  discipline  is  concerned 
with  choosing  the  geometry  or  composition  of  a  load-bearing  structure  so 
as  to  use  the  available  materials  as  efficiently  as  possible.  The  subject  has 
a  rich  history  and  an  extensive  literature;  books  and  articles  reviewing 
various  aspects  include  Refs.  7,  21,  53,  57.  Initially  attention  was  focused 
primarily  on  analytical  methods — optimality  conditions,  conformal  mapp¬ 
ing,  isoperimetric  inequalities,  and  so  forth.  More  recently,  with  the 
growing  feasibility  of  large  scale  computing,  attention  has  naturally  been 
turned  to  methods  for  the  direct,  numerical  calculation  of  optimal 
structures. 

To  fix  ideas,  let  us  consider  a  particular  problem  involving  shape 
optimization  and  plane  stress.  We  begin  with  a  homogeneous,  isotropic 
elastic  body  occupying  a  region  loaded  along  its  boundary  d(l  by 

a  specified  traction  /  We  desire  to  lighten  this  body  by  removing  material 
from  a  subset  //  c  fl,  consisting  of  one  or  more  holes  of  arbitrary  size  and 
shape.  The  goal  is  to  achieve  the  minimum  possible  weight,  i.e.  to  maximize 
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the  area  of  the  ‘holes’  H,  subject  to  a  performance  constraint  on  the  stress 
or  displacement  vi„  of  the  resulting  elastic  structure.  Typical  constraints  are 

that  the  work  done  by  the  load  (‘compliance’) 

be  not  too  large;  or  (5.1a) 

Jaa 

that  the  average  displacement  on  a  subdomain 

fij  be  not  too  large;  I  lu^l^Cor  (5.1b) 

Joi 

that  the  pointwise  maximum  stress  be  not  too 

large;  sup  ||  a„(x)  ||  ^  C  (5.1c) 

Highly  efficient  and  sophisticated  algorithms  have  been  developed  for 
the  numerical  solution  of  such  problems;  Ref.  21  gives  an  excellent  review. 
Typically,  one  begins  by  deciding  how  many  holes  to  consider.  Each  hole 
boundary  is  determined  by  finitely  many  points,  for  example  using  splines. 
The  resulting  domain  is  triangulated,  and  the  equations  of  elastostatics  are 
modeled  as  a  finite  system  of  linear  equations  using  the  finite  element 
method.  The  design  problem  is  thus  transformed  to  a  (highly  nonlinear!) 
mathematical  programming  problem,  and  one  can  seek  an  ‘optimal’ 
design— or  at  least  an  improvement  of  a  given  design — using  steepest 
descent,  or  perhaps  some  more  sophisticated  method. 

Though  its  utility  is  beyond  dispute,  this  ’conventional’  approach  has  one 
troublesome  aspect;  the  gross  features  of  the  design — especially,  the 
number  of  holes — must  be  chosen  at  the  outset;  they  are  not  a  part  of  the 
optimization.  Thus  the  output  is  likely  to  be  a  local  optimum,  or  at  best  an 
optimum  among  all  designs  with  a  specified  number  of  holes.  In  fact, 
numerical  attempts  at  global  optimization  for  related  model  problems  have 
led  in  some  cases  to  ‘optimal’  designs  that  vary  on  the  scale  of  the  mesh  size 
itself,  with  no  convergence  evident  as  the  mesh  size  tends  to  zero!*'  This 
phenomenon  is  now  well  understood.  In  the  context  of  shape  optimization, 
the  situation  is  as  follows;  consider  first  the  best  design  with  one  hole,  then 
that  with  two,  and  so  forth.  As  the  number  of  holes  gets  larger  the 
performance  may  get  better  (depending  of  course,  on  the  specific  problem 
under  consideration).  In  the  limit  of  infinitely  many  holes  one  thus  finds 
a  global  optimum  which  is  not  a  ‘conventional’  design  at  ail,  but  instead 
a  structure  made  from  composite  materials  obtained  by  perforation. 
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With  hindsight  it  seems  almost  obvious:  if  one  is  prepared  to  consider 
designs  with  many  small  holes,  then  one  ought  also  to  consider  their  limits. 
We  thus  arrive  at  a  new  approach  to  structural  optimization:  if  the  goal  is  to 
find  a  global  optimum  then  it  is  best  to  work  from  the  start  in  the  class  of  all 
structures  made  up  of  composite  materials  obtainable  by  perforation  from  the 
one  initially  given.  It  should  be  emphasized  that  the  underlying  problem  is 
not  being  changed,  since  we  allow  only  composites  achievable  by  perfora¬ 
tion,  and  we  are  careful  to  model  them  properly.  However,  the  resulting 
optimization  problem  looks  quite  different:  whereas  initially  we  were 
considering  structures  made  up  of  a  single  material  (or  the  absence  thereof), 
now  we  propose  to  allow  a  continuum  of  materials — each  representing 
a  perforated  composite  with  a  different  microscopic  geometry.  (As  a  tech¬ 
nical  matter,  the  mathematical  theory  discussed  in  the  preceding  sections 
does  not  quite  apply  to  perforated  composites,  since  it  requires  /ij>0  and 
Kj>0.  This  can  ^  circumvented,  at  least  for  compliance  optimization 
problems,  by  the  methods  of  Refs.  30, 32.  Alternatively,  we  can  simply  treat 
the  ‘holes’  as  though  they  were  filled  with  a  very  weak  elastic  material.) 

The  introduction  of  composites  as  generalized  designs — sometimes 
called  the  relaxation  of  the  design  problem — has  been  studied  extensively 
by  several  groups  over  the  past  ten  years,  see  for  example.  Refs.  20,  30,  38, 
50, 52, 54, 58, 59, 67.  From  a  theoretical  standpoint,  the  principal  advantage 
of  relaxation  is  that  it  assures  the  existence  of  an  optimal  design;  roughly, 
this  means  that  a  numerical  solution  of  the  relaxed  problem  will  converge  as 
the  mesh  size  tends  to  zero.  There  is  also  a  practical  advantage,  based  on  the 
fact  that  the  initial  material  and  the  absence  of  material  are  included  (as 
extreme  cases)  among  the  candidate  composites:  evidently,  for  a  given  finite 
element  subdivision  the  introduction  of  composites  serves  to  enlarge  the 
design  space  and  hence  to  improve  the  performance  of  a  numerically 
obtained  optimal  design.  Moreover,  precisely  because  it  has  the  effect 
(within  a  finite  element  context)  of  enlarging  the  design  space,  the  process  of 
relaxation  can  destroy  local  minima — making  it  easier  to  locate  a  globally 
optimal  design.  Finally,  since  the  relaxed  problem  is  known  to  have 
a  solution,  it  is  meaningful  to  use  the  associated  optimality  conditions;  this 
has  led  in  some  contexts  to  closed-form  examples  of  optimal  designs 
making  use  of  composites,  for  example.  Refs.  29-31.  The  method  of 
relaxation  has  its  limitations:  the  optimal  designs  obtained  this  way  may  be 
difficult  or  even  impossible  to  manufacture,  because  of  the  presence  of 
composites.  Even  so,  these  solutions  can  be  used  as  benchmarks  against 
which  to  compare  the  output  of  a  more  conventional  algorithm. 

The  process  of  relaxation  is  conceptually  simple:  we  must  simply 
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reformulate  the  design  problem  in  a  form  that  permits  perforated 
composites  as  admissible  materials.  The  actual  execution,  however,  is  not 
so  simple:  it  requires  specific  knowledge  about  the  properties  of  the  relevant 
composites.  For  a  local  performance  criterion  such  as  the  maximum  stress 
(S.lc)  we  would  have  to  know  optimal  bounds  relating  the  effective  Hooke’s 
law,  the  density  of  holes,  the  average  stress,  and  the  local  maximum  stress  in 
a  general  perforated  composite.  This  represents  a  challenge  for  the  future: 
no  such  result  is  presently  known.  For  a  performance  criterion  involving 
some  integral  of  the  displacement,  such  as  (5.16),  it  would  suffice  to  know 
the  solution  of  the  G-closure  problem — in  other  words,  to  know  the  class  of 
all  effective  Hooke’s  laws  obtainable  using  perforations  that  remove  a  given 
fraction  of  the  material.  The  analogous  problem  has  been  solved  for  scalar 
equations,^^  **  and  it  has  been  applied  to  solve  various  optimization 
problems  involving  conductivity,  see  for  example,  Refs.  1 2, 20, 30, 38, 50, 67; 
but  unfortunately  the  G-closure  problem  for  elasticity  remains  open  at  this 
time  except  in  certain  rather  special  cases.  However,  problems 
involving  compliance  constraints  such  as  (5.1a)  do  not  require  the  full 
solution  of  the  G-closure  problem;  rather,  bounds  of  the  type  presented  in 
Sections  3  and  4  are  sufficient.  To  explain  why,  we  note  that  it  is  not  really 
necessary  to  consider  all  composites;  one  might  as  well  consider  just  those 
that  can  actually  occur  in  an  optimal  design.  Now,  by  Green’s  formula  the 
compliance  is  equal  to  the  internal  elastic  energy: 


(5.2) 


{A(x)e{u),e{u))dx 


where  A{x)  is  the  spatially  varying  tensor  of  elastic  moduli  and  u  the 
associated  displacement.  A  structure  which  minimizes  weight  for  fixed 
compliance  will  also  minimize  compliance  for  given  weight;  it  is  not  hard  to 
see  from  this  that  A(x)  should  maximize  (Ae{u),e{y))  at  each  point  x  in  an 
optimal  design.  Thus  the  values  that  A{x)  can  take  in  an  optimal  design  are 
restricted  to  those  that  maximize  {A^,  for  some  tensor  cj. 

The  preceding  discussion  shows  that  we  have  enough  information  to 
solve  optimal  design  problems  with  compliance  constraints,  but  it  falls 
short  of  specifying  an  algorithm  to  do  so.  How,  operationally,  should  one 
proceed?  Following  Ref.  30,  we  advocate  an  algorithm  based  on  the 
principle  of  minimum  complementary  energy,  a  variational  principle  for  the 
stress  whose  extremal  value  is  equal  to  the  compliance: 


(5.3) 
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Introducing  a  Lagrange  multiplier  for  the  performance  constraint  (5.1a) 
our  design  problem  is 

(WEIGHT  +  A  •  COM  PLIANCE}  (5.4) 


The  outer  minimization  over  designs  is  quantified  by  introducing  functions 
6{x)  and  A{x),  the  density  and  effective  Hooke’s  law,  constrained  by  the 
pointwise  conditions 


0^0$  1,  and  A  is  the  effective  Hooke’s  law  of  a 

perforated  composite  obtained  by  removing  volume 

fraction  1  -  0  of  the  initially  given  material.  (5.5) 


The  compliance  is  itself  a  minimum,  according  to  (5.3),  so  (5.4)  becomes 


mm 

e,A 


0(x)dx  +  A*  mm 

O  div  a 


n  I  ( 


(A  ,  (x)CT,(T)dx 


(5.6) 


The  order  of  minimization  is  unimportant,  and  switching  it  gives 


with 


rnin  I  O^fff)  dx 

(5.7) 

div  a^0,o-n~f  JQ 

<I>^((T)  =  min  [0  +  x(A  ' '  cr,  cr)] 

e.A 

(5.8) 

The  minimization  in  (5.8)  is  over  real  numbers  6  and  tensors  A,  constrained 
by  (5.5).  This  is  slightly  different  from  the  problem  we  treated  in  Sections 
3  and  4,  but  it  can  be  solved  by  exactly  the  same  method — as  can 
considerably  more  general  problems,  for  example  the  analogue  of  (5.8) 
when  there  are  compliance  constraints  under  two  or  more  loads. 

The  next  step,  of  course,  is  to  evaluate  (5.8)  analytically  or  numerically, 
and  to  carry  out  the  optimization  by  solving  (5.7)  for  realistic  design 
problems.  Work  in  these  directions  is  currently  in  progress.  The  minimiza¬ 
tion  of  (5.8)  was  executed  in  Ref  30  for  the  special  case  of  an  elastic  material 
in  plane  stress  with  Poisson’s  ratio  zero,  i.e.  when  ix  =  k  =  \E,  where  E  is 
Young’s  modulus,  using  a  different  method,  based  on  quasiconvexification. 
The  answer  is  surprisingly  simple:  scaling  A  =  £  =  1  for  simplicity, 

l2(|(7i|-(-|ff2|)-2|CTi<72|,  |(T,|-|-|CT21^1 


where  ff,  and  are  the  principal  stresses  (the  eigenvalues  of  a). 
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